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1 Introduction 


A detailed account of progress achieved under this grant funding is contained in five journal 
papers. The titles and authors of these papers are: 

• Schols, J. and E. W. Eloranta, 1992: The calculation of area-averaged vertical profiles 
of the horizontal wind velocity using volume imaging lidar data, Journal of Geophysical 
Research, 97, 18395-18408 

• Eloranta, E. W. and D. K. Forrest, 1992: Volume imaging lidar observation of the 
convective structure surrounding the flight path of an instrumented aircraft, Journal 
of Geophysical Research, 97, 18383-18394. 

• Piironen, A. and E. W. Eloranta, 1995: Convective boundary layer mean depths, cloud 
base altitudes, cloud top altitudes, cloud coverages, and cloud shadows obtained from 
Volume Imaging Lidar data. Accepted for publication in the Journal of Geophysical 

Research. 

• Piironen, A. and E. W. Eloranta, 1995: An accuracy analysis of the wind profiles 
calculated from Volume Imaging Lidar data. Accepted for publication in the Journal 
of Geophysical Research. 

• Young, P. and E. W. Eloranta, 1995: Calculation of divergence and vertical motion 
from volume-imaging lidar data, Accepted for publication in the Journal of Geophysical 

Research. 

Copies of these papers form the body of this report. 
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Calculation of Area- Averaged Vertical Profiles of the Horizontal 
Wind Velocity From Volume-Imaging Lidar Data 

J. L. SCHOLS 1 AND E. W. ELORANTA 


Department of Meteorology, University of Wisconsin , Madison 

Area-averaged horizontal wind measurements are derived from the motion of spatial inhomogene- 
ities in aerosol backscattering observed with a volume-imaging lidar. Spatial averaging provides high 
precision, reducing sample variations of wind measurements well below the level of turbulent 
fluctuations, even under conditions of very light mean winds and strong convection or under the 
difficult conditions represented by roll convection. Wind velocities are measured using the two- 
dimensional spatial cross correlation computed between successive horizontal plane maps of aerosol 
backscattering, assembled from three-dimensional lidar scans. Prior to calculation of the correlation 
function, three crucial steps are performed: (1) the scans are corrected for image distortion by the wind 
during a finite scan time; (2) a temporal high pass median filtering is applied to eliminate structures that 
do not move with the wind; and (3) a histogram equalization is employed to reduce biases to the 
brightest features. 


1 . Introduction 

Convective motions in an unstable atmospheric boundary 
layer produce fluctuations in wind speed and direction. 
These velocity fluctuations scale with the convective scaling 
velocity (w*) and have typical magnitudes of 1 m s' 1 
[Kaimal et al ., 1976], which are spatially correlated within 
individual thermals. Thermal sizes scale on the boundary 
layer depth and have typical dimensions of the order of 1 km; 
when longitudinal rolls exist the stream-wise correlation 
length becomes much larger. These velocity fluctuations 
limit the accuracy with which the mean boundary layer wind 
can be estimated. To reduce errors in the mean much below 
w* , instruments that make point measurements must rely on 
temporal averaging over a period that is long compared to 
the lifetime of a thermal or long compared to the time for a 
thermal to advect over the instrument. When the mean wind 
is small, the required averaging times become long. For 
example, with alms -1 mean wind and a 1-km-deep 
boundary layer, a single cell requires approximately 1000 s 
to pass a fixed sensor. Reduction of the sampling error to a 
level much below w* requires the passage of many thermals; 
this in turn requires many hours of averaging. 

Spatial averaging offers the possibility of reducing ther- 
mal-scale fluctuations in the wind without long averaging 
times. This paper describes a lidar technique for obtaining 
area-averaged winds measured using lidar observations of 
the drift of inhomogeneities in naturally occurring atmo- 
spheric aerosols. Spatial inhomogeneities in atmospheric 
boundary layer aerosol content are produced by a variety of 
processes, including localized sources, and by the vertical 
entrainment of air from above. The displacement of these 
inhomogeneities can be measured between successive ob- 
servations by an incoherent lidar system. These displace- 
ments allow calculation of wind velocities without the use of 
Doppler techniques. Wind measurements have been ob- 
tained by calculating correlations between a time sequence 
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of lidar profiles acquired from a small number of discrete 
lidar printing directions [Eloranta et al. , 1975; Kunkel et al ., 
1980; Sroga et al . , 1980; Hooper and Eloranta , 1986; Kolev 
et al., 1988] and by calculating cross correlations from a 
sequence of two dimensional lidar images [Leese and Ep- 
stein , 1963; Sasano et al ., 1982; Ferrare et al., 1991]. 

The University of Wisconsin (UW) volume imaging lidar 
(VIL) is designed to produce a continuous sequence of 
three-dimensional maps of aerosol backscattering [Eloranta, 
1987, 1988]. A relatively fast laser pulse repetition rate, fast 
angular beam scanning, and a fast data acquisition system 
allow this lidar to image atmospheric boundary layer (ABL) 
structure over an approximately 50- to 70-km 2 area at 
approximately 3- to 4-min intervals. Each volume image 
contains approximately 6 million measurements of aerosol 
backscattering. 

In the present paper, lidar profiles observed with the VIL 
are used to produce images of the aerosol structure on 
vertically stacked horizontal planes. These images are com- 
pared to similar images obtained a short time later to obtain 
independent estimates of the wind velocity at the altitude of 
each plane. The vertical profiles of the horizontal wind 
obtained by this process represent an area average over ~ 50 
to 70 km 2 . 


2. Method of Data Analysis 

The backscattered power incident on the lidar receiver 
from a pulse of laser light is described by the single- scatter 
lidar equation [Collis and Russell, 1976], which, in a slightly 
modified form, reads 


E A c 

P(r) = -j ~ P'(n, r) exp 
r 2 


where 



r = range from lidar; 

P(r) = received power; 

E = transmitted laser pulse energy; 

A - area of receiver; 
c- speed of light; 

/3'(7r, r) = volume backscatter coefficient per steradian, 
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which is equal to the product of the aerosol 
scattering cross section, per unit volume and 
the aerosol backscatter phase function; 
ft e (r) = extinction cross section per unit volume. 

Before recording lidar profiles on the optical disk, the VIL 
data system (1) normalizes each lidar profile to remove the 
energy variations between laser pulses, (2) corrects for the 
inverse range squared dependence; and (3) forms the natural 
logarithm of the normalized lidar profile. The recorded 
logarithm of the range corrected energy normalized data 
points Sy are given by 

(2) 

where 

rj = 8rj\ the distance between the lidar and yth scattering 
range element; 
i= index for tth lidar profile; 
j = index for yth range element; 

5r = range increment between successive data points. 

As a first step in the data analysis each lidar profile is 
high-pass filtered by subtracting a moving median of data 
points in a range-centered interval around the jth range 
element. The filtered lidar profile, s{j , is given by 

S{j=S 0 -S^ (3) 

where 

S™ = median {S i} } j - n < / <j + n (4) 

and In is the width of the spatial median. This filter 
eliminates small intensity fluctuations produced by errors in 
the energy normalization of individual lidar profiles and 
large-scale spatial variations produced by atmospheric atten- 
uation. A median filter is used instead of a running average 
because it is much less sensitive to perturbation by extreme 
values, such as cloud reflections. In the ABL, aerosol 
inhomogeneities, which have spatial scales of the order of 
the boundary layer depth typically have the largest backscat- 
ter fluctuation [Kunkel et al ., 1980]. Since the boundary 
layer depth during the present experiment was approxi- 
mately 1 km, a filter length of 2 km proved to be a suitable 
choice. 

The filtered lidar profiles making up a volume scan provide 
a spherical coordinate (range, elevation, azimuth) represen- 
tation of aerosol inhomogeneities. In the second step of the 
data analysis the representation is changed to cylindrical 
coordinates (altitude, range, azimuth). In this new constant 
altitude plan position indicator (CAPPI) format, individual 
profiles appear to have been acquired by a lidar located at an 
altitude, z, and scanning horizontally through the same 
azimuth angles sampled in the original data. These new 
profiles are formed by linear interpolation between the 
nearest two lidar profiles recorded at the same azimuth 
angle. Figures 1 and 2 show the geometry of this conversion. 

Wind velocities are obtained from the CAPPI scans using 
an algorithm that includes the steps described below. 

2.1. Correction for Wind Distortion 

During the 3 to 4 min required to complete a volume scan, 
the aerosol structures move along with the wind. The spatial 


structure is therefore distorted in the resulting lidar images. 
This effect is illustrated in Figure 2. During one scan the 
scanning radius, r, completes an arc, C, of length r(d<£/ 
dt)T s , where (d<t>/dt) is the azimuth angular scanning rate, <f) 
is the azimuth angle, t is time, and T s is the time needed to 
complete a volume scan. At the same time, an aerosol 
structure moves over a distance D - V t T s , measured along 
the scanning direction (V, is the component of the aerosol 
velocity perpendicular to the scanning radius). 

The shape distortion of the aerosol structures in the 
CAPPI scans is corrected by a method similar to that 
described by Sasano et al. [1982]. The correction for shape 
distortion is carried out by moving the individual aerosol 
backscatter profiles making up a CAPPI scan a distance 

d = v • 6/, (5) 

where v is the average horizontal wind velocity in the scan 
area and 8t is the time difference between the start of a 
CAPPI scan and the time when this azimuth is sampled. 
Because profiles making up the CAPPI scan have been 
computed from a number of profiles in an elevated scan, 8t 
is approximated with the time difference from the start of the 
volume scan to the start of a range-height indicator (RHI) 
scan at this azimuth. 

Finally, the data in each CAPPI scan were transformed 
from a cylindrical coordinate system to an image on a 
Cartesian x, y grid, in which the x axis was taken along the 
middle azimuth direction of a CAPPI scan. The origin of the 
x, y grid coincides with the VIL position. The data values on 
the Cartesian grid points were found through weighted linear 
interpolation between the grid points in the cylindrical 
coordinate system. The resolution of the Cartesian grid, 
typically approximately 45 m, was chosen to appropriately 
match the spatially averaged resolution of the cylindrical 
coordinate system. 


2.2. Temporal High Pass Median Filtering 

The lidar images may contain structures that do not move 
with the wind. Such structures typically represent residual 
effects of attenuation, small nonlinearities in the lidar system 
response, or atmospheric features, which are anchored to 
the underlying terrain. To remove these features, a temporal 
median image, / m (x, /,), is constructed at each level, over a 
time interval centered about each image 

f m (x, tf) = median {/(x, tj)} (6) 

/ - 1 < 7 < i + 1 , 

where /(x, tf is the instantaneous aerosol back scattering 
distribution in an image at time /,• and at grid location x = (x, 
y), and 2/ is the width of the temporal median. The median 
image is thereafter subtracted from each image to obtain the 
temporal high pass median filtered image, /^(x, tf: 

/(x, tj) =f(\, tj) tj). (7) 

Ideally, the median image is computed over a time period 
that is short compared to the rate at which the stationary 
background image changes, yet long enough that many wind 
driven aerosol structures pass over each image point. In 
practice, filter durations between 30 and 60 min are typically 
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Fig. 1. Three of the 60 range-height indicator (RHI) scans, 
forming a complete volume scan, are shown. The scanning pattern 
follows the lines a-b-c-d-e. Successive RHI scans are separated by 
8<j> in azimuth. Synthetic constant altitude lidar profiles are then 
formed at selected altitudes, for example, along the dashed line h. 
Aerosol backscatter values on a constant-altitude line, between 1 
and 2, on a constant-altitude line, are found through linear interpo- 
lation between nearest data points in the lidar profiles p and q. 


employed. This process removes the stationary features 
without affecting the motion of the inhomogeneities. 

2.3. Histogram Equalization 

The instantaneous aerosol backscattering distribution, 
/(X, t), in an image can be represented as a distribution of 
pixel brightness values. Here x is the pixel position in x, y 
coordinates. Anomalously bright targets may dominate and 
thus bias the measured average velocity to the velocity of the 
brightest targets. To reduce this effect, the pixel brightness 
distribution is made uniform by changing pixel brightness 
values according to the transformation 

/„(x, t) = P[f f (x, ?)], (8) 

where / w (x, t) corresponds to the new pixel brightness 
distribution and P stands for the cumulative probability 
distribution of the initial pixel brightness distribution. This 
transformation is illustrated in Figure 3. Pixel brightness 
values that lie in the tails of the distribution are compressed, 
whereas the pixel brightness values in the central region are 
expanded. This operation is referred to as histogram equal- 
ization. 


Y 



Fig. 2. A CAPPI scan is corrected for wind distortion by 
displacing each of its radial aerosol backscatter profiles over a 
distance ~\8t, where 8t is the time elapsed since the start of a 
CAPPI scan. 



NORMALIZED BRIGHTNESS 


Fig. 3. (a) The probability density distribution (PDF), shown by 

the short-dashed curve, and the cumulative probability distribution 
(CPD), shown by the solid curve, of pixel brightness values for an 
image on July 1, 1987, at 1356:49 CDT at a height 400 m above the 
lidar site. ( b ) Same as in Figure 3 a, except after the histogram 
equalization, using 256 levels of pixel brightness. The transforma- 
tion for two pixel brightness values is illustrated in Figure 3 a with 
the dashed curves labeled a and b. 


2.4. Spatial Cross Correlation Function 

To find an objective estimate of the aerosol motion aver- 
aged over the area of an image, the location of the maximum 
value of the two-dimensional spatial cross correlation func- 
tion (CCF) is used. The CCF between the patterns /„(x, t x ) 
and / M (x, t 2 ), observed at the times f, and t 2 , respectively, 
is defined as 


C(Sx, t u t 2 - f|) = J / M (x, *j)/ M (x + 8x, t 2 ) dS, (9) 

where 8\ = (fix, fiy) is a lag vector and the integral is taken 
over an area S . 

The CCF is computed between successive images using 
the computationally efficient fast Fourier transform (FFT) 
technique [Otnes and Enochson , 1972]. Before the FFT is 
applied, each image is linearly detrended to take the spectral 
power out at zero wave number and to minimize leakage 
from low wave numbers into the higher wave numbers. The 
discrete FFT of an image on a regular x, y grid with 
dimensions L x and L y is written as 
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F„(k, 0=2 2 /«(*> ') exp' 0 ' 2 ”' 1 ’ 0 , (10) 

m = 1 n = 1 

where 

k = (A: Jt , k y ), the wave vector; 
k x = {i x M x )l(ML x ), the wave number in the x 
direction; 

i x =lag index in the jc direction {-MU ^ i x < Ml 2); 
M* = number of data points in the jc direction; 

M = number of data points in the jc direction, 
extended with zeroes to obtain a power of 2; 
k y ={i y N y )l(NLy), the wave number in the y 
direction; 

i y =lag index in the y direction (-N/2 < i y < Nl 2); 
N y = number of data points in the y direction; 

N= number of data points in the y direction, 
extended with zeroes to obtain a power of 2; 

j - v=h 

X = (jc, y), the coordinate of a data point on a jc, y 
grid or ( mL x IM x , nL y lN y )\ 

L x fM x =grid resolution in the jc direction; 

L y /Ny =grid resolution in the y direction. 

The basic method uses the FFTs of two successive images 
to compute the cross spectrum, G(k, t u t 2 ~ t j), where 

G(k, t u ti~ t\) = 4 ~rp> F*(k, /,)F M ( k, t 2 ) (11) 

M x Ny 

and an asterisk denotes a complex conjugate. The cross 
spectrum is then high-pass filtered in wave number space to 
eliminate large spatial scale coherence, which may be caused 
by residuals introduced by the temporal high pass filtering 
process. Care must be taken here that the spatial high pass 
filtering process does not remove the large-scale coherent 
aerosol structure between successive images. For a typical 
ABL, where convective structures scale on the boundary 
layer depth, ~1 km, a filter that removes scales larger than 2 
km by 2 km is a good choice. The inverse FFT of the cross 
spectrum results in the cross-covariance function, 

1 M x Ny 

Cov (5x, = - — — 

M N 

• 2 2 '2 - h) exp 0*2 TTk • 8x). (12) 

/ = l / = l 

The covariance function is the numerator for the definition of 
the CCF, 


where 

y m x Ny 

f«.ave(*, *) = 2 2 /«(*’')• ( 15 > 

x y m = 1 n = 1 

To reduce statistical error, an average correlation function, 
C ave (Sx, tj, tj - ti), is then computed over successive 
images 

1 3 ~ 1 

tj , tj t /) — ” 7 C{8x, t( + | ~ t () 

(16) 

where the indices i and j refer to the times between which 
the averaging is carried out. 

2.5. Location of the Correlation Maximum 

To estimate the wind speed, it is necessary to locate the 
maximum of the average correlation function, C ave (Sx, tj , tj 
- tj). Interpolation is required because this function is only 
available at discrete points, 8x = (mL x /M xy nL y IN v ). In 
addition, the correlation function is also often quite flat near 
its peak, so small noise fluctuations can move the apparent 
peak to one of several points adjacent to the maximum. An 
approximate location for the maximum, 5x max , is first ob- 
tained by searching for the peak in the two-dimensional 
array C ave (Sx, t Jy tj - tj) y which is smoothed by using a 
moving area average over 3 by 3 points. A bicubic natural 
nonsmoothing spline interpolation is applied to the un- 
smoothed average correlation function around the array 
location <5x max . The position of the maximum of the spline 
function determines the horizontal displacement of the aero- 
sol inhomogeneities. Dividing this displacement by the time 
separation between two successive images gives the wind 
velocity. A bicubic natural nonsmoothing spline function is 
used for the interpolation around 6 x max , since that function is 
smooth and preserves the shape of C ave (5x, t jy tj - tj). Use 
of the bicubic natural nonsmoothing spline interpolation 
improves the accuracy of the determination of the position of 
the maximum of the average correlation function to within a 
fraction of the grid resolution. 

A first estimate of the horizontal wind velocity is obtained 
using this algorithm without the correction step for wind 
distortion. The first estimate of the wind velocity is then 
used to repeat the algorithm in its full extent to get a better 
value of the wind velocity. The algorithm converges rapidly, 
requiring only one or two iterations. 

3. Results 


C(5x, r,. 


*2-'l) = 


Cov (8x, t j , t 2 - U) 

\ cr ] (j 2 


(13) 


where cr, and a 2 are the standard deviations of the succes- 
sive images. The square of the standard deviation of an 
image is defined as 


cz 2 = 


1 

{M x Ny - 1 ) 


Mx Ny 


2 2 -/«. ave(x> t )] 2 , 


m = 1 n = 


(14) 


3.1. Site and Instrumentation 

Measurements were obtained between June 25 and July 
15, 1987, on the Konza Prairie long-term ecological research 
(LTER) site near Manhattan, Kansas, as part of the First 
International Satellite Land Surface Climatology Project 
(ISLSCP) Field Experiment (FIFE). Sellers et al. [1988] give 
a brief description of the experiment. The primary 15-km by 
15-km field site is a mostly treeless, rolling prairie with 
altitude variations of roughly 100 m. The UW VIL system 
was located near to the center of the site (96°32'30"W, 
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Fig. 4a. Perspective three-dimensional view of convective cells in the atmospheric boundary layer (ABL) 
displayed over a high-resolution SPOT picture of surface topography on July 1, 1987, at 1121 CDT under almost 
windless convective conditions. This display is prepared by selecting a contour layer to construct a solid surface 
through the volume of normalized aerosol backscatter data. The lidar is located at the origin of the coordinate system 
provided in the picture. All dimensions are given in kilometers. A fork-shaped creek bed area is clearly visible on the 
surface, with one leg at 8 km and the other one between 5 and 8 km. 


TABLE 1. University of Wisconsin Lidar System 
Parameters (1987) 


Instrument 

Parameter Specification 

Laser transmitter 

wavelength 

1064 nm 

beam pointing resolution 

0.2 mrad 

average output energy per shot 

0.6 J 

pulse repetition rate 

30 Hz 

pulse duration 

~7 ns 

Receiver telescope 

diameter 

0.5 m 

field of view 

adjustable, 1.5-7 mrad 

detector 

Silcon Avalanche 
Photodiode 
**(RCA 30955E) 

quantum efficiency 

40% at 1064 nm 

spectral band pass 

1.0 nm, interference filter 

Data logging 

logarithmic amplifier 

100-)nV to 6-V input range 
(=90dB) 

A/D resolution 

10-bit resolution 

sampling rate 

20 MHz max, 10 MHz this 
study 

range resolution 

7.5 m, 15 m this study 

maximum range 

15 km, this study 

preprocessing 

VAX 11/750, LSI 11/73, and 
CSPI array processor 

output 

2.6 GByte write once optical 
disk, real-time display 


39°03'51"N) at 448 m above mean sea level (msl). The 
average elevation of the FIFE site is approximately 408 msl. 
A summary of the UW three-dimensional VIL system pa- 
rameters is presented in Table 1. A Doppler lidar system, 
operated by NOAA, was located next to the UW VIL. 
Sixteen portable automated mesonet II (PAM II) stations 
operated by NCAR were located within 10 km of the UW 
VIL. Members of Cornell University also operated a radio- 
sonde system. The balloon launch location was at 
96°33'47"W, 39°07'02"N, at 340 msl, approximately 6 km 
north-northwest of the lidars. 

3.2. Description of the Data Acquisition 

The VIL provided a time sequence of three-dimensional 
maps of aerosol backscattering. Typical lidar volume scans 
consisted of 60 RHI scans between the horizon and 10° of 
elevation, with about 60 lidar profiles, each consisting of 
1024 range gates, in each RHI. Individual RHI scans were 
spaced 0.5° apart in azimuth to cover a 30° azimuthal sector. 
This resulted in a time separation of 165 s, between two 
successive images, which was short compared to the lifetime 
of the coherent aerosol structures. It was necessary to 
compromise between the angular width of the volume scans 
and the frequency of observation of the scans, as con- 
strained by the laser repetition rate. The angular size and 
timing of the images limits the maximum wind speed that can 


A/D is analog to digital. 
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Fig. 4b. CAPPI scan constructed 500 m above the lidar site using the data on July 1. The bright areas correspond 
to aerosol-laden structures that originate from the surface, whereas the dark regions indicate clear air between these 
structures. 


be measured. If the scan time is too long or the angular 
sector is too small, the aerosol structure observed on one 
scan may have advected out of the scan area by the start of 
the next scan. 

Two data sets, acquired under very different conditions, 
were chosen to test the wind-measuring algorithm. The first 
set was collected between 1321 and 1424 CDT on July 1 
(CDT is UT - 5 hours) in a convectively driven boundary 
layer. The second data set was acquired between 0627 and 
0657 CDT on July 7 in a largely shear driven boundary layer. 
Representative volume and CAPPI scans are shown in 
Figures 4 and 5. Both data sets show bands of organized 
aerosol structures. On July 1 these bands were mainly 
aligned with the wind. The weather was clear, with scattered 
convective clouds. Surface winds measured by PAM sta- 
tions were less than 2ms" 1 . The depth of the ABL changed 
from approximately 1 100 m above the lidar site at the start of 
the data set to around 1300 m at the end. Scattered convec- 


tive clouds appeared on top of the convective plumes. On 
July 7 the aerosol bands were orientated nearly perpendic- 
ular to the wind direction. The sky above the site was clear, 
and the surface sensible heat flux was small. A thunderstorm 
was approaching, and the winds were 3 to 10 m s -1 . This 
data set provided an example of shear-driven turbulence, in 
contrast to the convective conditions of July 1. Boundary 
layer depth grew from 300 m at the beginning of the period to 
500 m at the end. The contrast between the aerosol struc- 
tures and their background, that is, the image contrast, 
decreased significantly during this period. 

3.3. Computational Results and Accuracy 

Figures 6a-6d illustrate the effect of various processing 
steps applied to the CAPPI scans obtained on July 1. The 
convective, low wind speed conditions provided a particu- 
larly good test case for the performance of the temporal high 




Same as in Figures 4 a and 4 b, except for a shear-driven ABL with low surface heat flux on July 7, 1987, at 
0636 CDT. The CAPPI scan is made at 250 m above the lidar site. 
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4 6 8 10 12 14 

Distance, km 

Fig. 6. The temporal high pass median filtering step and the histogram equalization step for the July 1 data set at 
1321:09 CDT at 100 m above the lidar. (a) CAPPI scan. ( b ) The temporal median, taken over 24 images (66 min), 
displays the motionless background. The fork-shaped dark area around 8 km reflects the creek beds on the surface, (c) 
The temporal high pass median filtering removes the motionless background. ( d ) The image after histogram 
equalization. 


pass median filtering step. The median image shows a clear 
picture of the background in which the aerosol structures are 
embedded. The histogram equalization step, which was 
applied to make the brightness distribution in the image 
uniform, enhances the visual contrast in the gray regions. 


Figure 7 displays a comparison of VIL-determined winds 
with conventional measurements. On July 7 all profiles (top 
panels of Figure 7) show similar vertical shear. The UW VIL 
profile at 0627:18 CDT and the Doppler lidar profile at 0611 
CDT agree very well, except below 300 m height. The UW 
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Fig. 7. A comparison of wind profiles measured by the UW VIL, the NOAA Doppler lidar, and the radiosonde. 
(Top panels) Wind profiles on July 7, 1987, where the solid curve shows the UW VIL wind profile at 0627 CDT; the 
short-dashed curve shows the Doppler lidar wind profile at 0611 CDT; the long-dashed curve shows UW VIL wind 
profile at 0654 CDT; and the dotted-dashed curve shows the radiosonde wind profile at 0704 CDT. The UW VIL wind 
profile is calculated at altitudes of 15, 30, 60, 100 m and then at 50-m intervals above 100 m. The radiosonde profile is 
averaged over 100-m layers. (Bottom panels) Wind profiles for July 1, 1987, where solid curve shows the 660-s average 
UW VIL wind profile at 1326 CDT; the short-dashed curve shows the Doppler lidar wind profile at 1327 CDT; the 
long-dashed curve shows the 660 s average UW VIL wind profile at 1337 CDT; and the dotted-dashed curve shows the 
radiosonde wind profile at 1339 CDT. The UW VIL wind measurements were calculated at 100-m vertical intervals. 
Thirty-minute averaged 10-m altitude winds measured by portable automated mesonet (PAM) stations are also shown. 
Two separate 30-min averages are shown for each PAM station; the first average (circles) begins at 1300 CDT. The 
second average period (crosses) begins at 1330 CDT. 


VIL profile at 0654:44 CDT and the vertically averaged 
radiosonde profile at 0704 CDT, match closely. The UW VIL 
profile looks smoother, as a result of the averaging over a 
32-km 2 horizontal area. In the early afternoon of July 1 
(bottom panels of Figure 7) the conditions were convective, 
with almost no average wind. These conditions, together 
with the organized roll convection (see Figure 4), caused 
large updrafts and downdrafts in the ABL. As a result, bright 
and dark bands appeared in the images. To obtain a stable 
estimate of the horizontal wind velocity, it is therefore 
necessary to average over an area of 12 km by 6 km. This 
allowed averaging over many convective cells to obtain a 
consistent estimate of the mean wind velocity. The fluctua- 
tions in a UW lidar profile between adjacent points are of the 
order of 0.1 ms -1 in speed and 10° in direction. As expected, 
under almost windless conditions with strong clear air con- 
vection, little vertical shear is evident. The Doppler lidar 
showed a wind speed that decreases with altitude. The 
difference between the UW VIL derived winds and the 
Doppler lidar determined winds are most likely caused by 
sampling effects, caused by the much smaller averaging 
volume represented by the Doppler velocity azimuth display 
(VAD) scan. The Doppler VAD wind measurement tech- 
nique provided wind measurements during a 2-min sweep for 
a volume consisting of a vertically oriented cone. The VIL 
was located at the vertex of the cone. Wind measurements at 


each altitude sample a ring around the lidar with a diameter 
that increases with altitude. The wind profile measured by 
the radiosonde balloon represents a point measurement of 
the wind velocity at each level. The wind profile measured 
by the radiosonde drift shows large variations that are 
caused by the turbulent motions on the balloon path. In this 
convective light wind case, the turbulent velocities scale 
with the convective scaling velocity w*, which was approx- 
imately 1.5ms' 1 . The 30-min average 10-m winds measured 
by PAM stations are also shown. Two successive averages 
are shown for each PAM station. Large differences exist 
between stations and between successive measurements. 
The scatter in the points on each profile, the differences 
between the profiles, and the variations in the PAM mea- 
surements are equal to or less than w*. 

Each of the individual steps of the wind measurement 
algorithm was found necessary to obtain smooth and accu- 
rate profiles. Figure 8 shows 165-s average wind profiles, 
averaged over an 8-km by 4-km area. They were acquired 
during the morning of July 7, between 0627:18 and 0630:03 
CDT, when the image contrast was the highest (top panels of 
Figure 8), and at the end of this data set, between 0654:44 
and 0657:29 CDT, when the image contrast was the lowest 
(bottom panels of Figure 8). The calculations were done 
using all steps and deleting one step each of the algorithm. 
The data set contained 12 images, and the time separation 
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Fig. 8. (Top) July 7, 1987, 165-s average wind profiles at 0627 CDT and (bottom) at 0654 CDT measured using all 
the steps in the algorithm (solid curve). The short-dashed curves contain no correction for wind distortion. The 
long-dashed curves have no temporal high pass median filtering; and the dotted-dashed curves have no histogram 
equalization. 
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Fig. 9. University of Wisconsin VIL wind profiles obtained on July 7 at 0654 CDT: (Top panels) the solid curve 
shows median length of 330 s and a temporal average of 165 s; the dotted curve shows a median length of 990 s and a 
temporal average of 165 s; the dashed curve shows a median length of 330 s and a temporal average of 330 s; and the 
dotted-dashed curve shows a median length of 990 s and a temporal average of 330 s. (Bottom panels) The 165-s 
averages calculated with and without spatial high pass filtering (see section 2.4) that eliminates scales larger than 2.5 km 
in both horizontal dimensions, where the solid curve shows the median length of 330 s with filtering and the dashed 
curve median length of 990 s with filtering; and the dotted-dashed curve shows the median length of 990 s and no spatial 
high pass filtering. 
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Fig. 10. Profiles measured in the afternoon of July 1, 1987, at 1326 CDT using two different lengths, namely, 2640 
and 3795 s, of the temporal high pass median filter and temporal averaging over 165 and 660 s. The short-dashed curve 
shows a median length of 2640 s and a temporal average of 165 s; the solid curve shows the median length of 3795 s and 
temporal average of 165 s; the long-dashed curve shows the median length of 2640 s and temporal average of 660 s; and 
the dotted-dashed curve shows the median length of 3795 s and temporal average of 660 s. 


between two successive images was 165 s. The correction 
for wind distortion has a larger effect in the upper part of the 
boundary layer, where the wind speed is higher. Wind 
direction is almost uniform throughout the boundary layer 
depth. Deleting the step at 0627:18 CDT leads to an overes- 
timation of the horizontal wind speed throughout the bound- 
ary layer, because the direction of the wind component 
transverse to the scanning radius coincides with the clock- 
wise scanning direction. At 0654:44 CDT, when there is little 
contrast between the aerosol structures and their back- 
ground, omission of the correction for wind distortion does 
not always lead to an overestimation of the horizontal wind 
speed. In spite of the histogram equalization step, a few local 
bright spots in a low contrast image may still dominate the 
CCF computations. The area-averaged wind velocity is then 
biased toward the apparent motion of these individual spots. 
Because the median background in which the aerosol struc- 
tures were embedded was changing rapidly, the length of the 
temporal high pass median filter was selected to be 330 s, 
twice the time separation between two successive images. 
The temporal high pass median filter most strongly affects 
the regions near the surface and at the top of the boundary 
layer. This is because topographic features are most promi- 
nent near the surface, while nonlinearities in the lidar system 
response or residual effects of attenuation are more likely 
near the top. When this step is not used, the wind speeds are 
biased to low values by stationary features. 

The computations showed little sensitivity to the length of 
the temporal high-pass median filter. Figure 9 shows the 
results, for the early July 7 data, after the image contrast had 
decreased substantially at the end of this data set. Different 
choices of the length of the temporal high pass median filter, 
namely, two and six time separations between two succes- 
sive images, resulted in differences smaller than 0.3 ms" 1 in 
wind speed and 10° in wind direction for a 165-s average 
profile (see top panels of Figure 9). These differences were 
reduced when more CCFs were averaged together. For a 
330-s average the differences between the profiles calculated 
for different temporal median lengths become smaller than 
0.1 m s _l in speed and less than 10° in direction. However, 
when the spatial high pass filtering of the cross spectrum (see 
section 2.4), which eliminates scales larger than 2.5 km, is 
omitted, the profile for a median length of six time separa- 


tions exhibits large errors (bottom panels of Figure 9). The 
high-pass spatial filter apparently reduces sensitivity to the 
length of the temporal median. Figure 10 shows the results of 
a calculation for different lengths of the temporal median, 
namely, 2640 and 3795 s, on July 1. For 660-s averages the 
difference between the profiles is smaller than 0.1 m s _1 in 
speed and 10° in direction. In addition, the profiles for the 
longer 24-image median look smoother, which indicates that 
this is a better choice for the median length for this convec- 
tive, low wind speed case. 

In Figure 11, CCFs are presented for both July 1 and July 
7 data sets. The CCFs are shown as three-dimensional 
surfaces to enhance their visualization. The height of the 
peak with respect to the background is clearly depicted. 


4. Discussion and Conclusions 

Horizontal wind speeds and directions within the ABL 
were estimated by measuring the drift of aerosol backscat- 
tering patterns inside successive lidar images, using a two- 
dimensional spatial cross correlation technique. 

The wind-measurement algorithm, described in section 2, 
has been validated during FIFE, under two very different 
conditions, namely, a convective light wind case and a 
shear-driven case with moderate to high wind speeds. The 
derived wind profiles were shown to be nearly insensitive to 
the choice of the length of the temporal high pass median 
filter, provided an additional spatial high pass filtering is 
used. The area averaging allowed accurate mean horizontal 
wind measurements in the presence of turbulence, which 
degrades conventional observations. 

A comparison of the lidar wind measurements with winds 
measured by Doppler lidar and balloon-borne radiosondes 
displayed differences of the order expected from the varia- 
tions produced by the different spatial and temporal averag- 
ing inherent in each technique. For example, on July 1, 
internal consistency between winds independently calcu- 
lated for adjacent altitudes or times suggest that spatial 
averaging reduced sampling errors well below the level of 
turbulent fluctuations. The variation between adjacent lev- 
els, of 5- to 10-min average lidar wind measurements, is 
smaller than 0.1 m s” 1 in speed and 10° in direction. 
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Fig 1 1 Cross-correlation functions (CCFs) displayed as three-dimensional surfaces for the July 1 and July 7 data 
sets. The axes are annotated in kilometers. (Top) July 1, 1321:09 to 1323:54 CDT at 100 m height above the lidar site. 
(Bottom) July 7, 0627:18 to 0630:03 CDT at 30 m height above the lidar site. 
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Volume-Imaging Lidar Observations of the Convective Structure 
Surrounding the Flight Path of a Flux-Measuring Aircraft 

Edwin W. Eloranta and Daniel K. Forrest 


Department of Meteorology , University of Wisconsin, Madison 

The University of Wisconsin volume imaging lidar has been used to portray images of the 
three-dimensional structure of clear air convective plumes in the atmosphere surrounding the flight 
path of the instrumented Twin Otter aircraft operated by the National Aeronautical Establishment 
(NAE) of Canada. Lidar images provide a context for interpretation of the aircraft measurements. The 
position of data points within a convective element can be determined and the temporal development 
of the plume can be observed to time the observation with respect to the life cycle of the plume. Plots 
of the vertical flux of water vapor, q'w ', superimposed on lidar images clearly demonstrate the 
well-known sampling difficulties encountered when attempting to measure fluxes near the top of the 
convective layer. When Loran was used to determine average aircraft velocity, flight-leg-averaged 
horizontal winds measured by the aircraft and area-averaged winds measured by lidar agree to within 
0,2 ms" 1 in speed and 1° in direction. 


Background 

Much of our knowledge of the planetary boundary layer 
has been acquired from aircraft-borne instrumentation. 
These observations have produced a rather complete de- 
scription of the mean properties of the mature, midafternoon 
convective boundary layer over homogeneous surfaces. 
Mean profiles of state variables as well as flux profiles for 
heat, momentum, and water vapor have been measured 
during many experiments and have allowed the development 
of a system of convective layer scaling that integrates this 
information into a usable form [Stull, 1988]. 

While aircraft instruments are well suited for observations 
of established boundary layers over homogeneous terrain, 
difficulties arise when the layer undergoes rapid temporal 
changes or the surface properties change rapidly in time or 
space. Convective layer scaling shows that the depth of the 
layer is a fundamental length scale and should be measured 
simultaneously with flux measurements inside the boundary 
layer. Interpolation errors are encountered when a single 
aircraft sequentially measures mixed layer depth and vari- 
ables inside the layer. To obtain statistically reliable esti- 
mates of boundary layer variables, aircraft measurements 
must be averaged over many turbulent eddies; this is partic- 
ularly true outside the atmospheric surface layer, where 
eddies have scales determined by boundary layer depth and 
flight legs with lengths equal to or exceeding tens of kilome- 
ters are often required [Lenschow and Stankov , 1986]. It 
clearly becomes difficult to separate effects of spatial vari- 
ability from statistical sampling variability. 

Even greater difficulties are encountered when aircraft 
data are used in attempting to describe the typical structure 
of individual turbulent eddies. Although precise in situ 
measurements can be obtained along the flight path, it is 
difficult to determine the aircraft position with respect to the 
eddy structure. Many efforts to explain eddy structure have 
relied on conditional sampling techniques, where an indica- 
tor function constructed from one or more of the measured 
variables is used to control the sampling of other variables in 
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a phase-coherent manner to generate a composite picture of 
a typical eddy [ Greenhut and Khalsa , 1987; Khalsa and 
Greenhut , 1987]. The picture of the eddy that results de- 
pends critically on the choice of the indicator function 
[Marht and Frank , 1988]. The measurement signature de- 
pends on how and when the aircraft penetrated the eddy. 
The measurement path may have penetrated through the 
center of the eddy or just grazed the edge. In addition, the 
timing of the observation with respect to the growth and 
decay of the eddy is unknown. It is very difficult to deter- 
mine whether a particular observation represents an eddy 
that is actively extracting energy from the large-scale shear/ 
buoyancy field or whether it is the inertial remains of a 
dissipating eddy. 

Remote sensors such as radar and lidar allow imaging of 
large volumes of the atmosphere in a short period of time, 
thus permitting visualization of atmospheric eddy structure 
[Hardy and Otterster , 1969; Kunkel et al. , 1977]. However, 
the images produced often lack the quantitative detail pro- 
vided by in situ sensing. Measurements normally consist of 
one or two variables. For example, radar or lidar reflectivity 
may be provided; if the system has Doppler capability, radial 
velocities may also be included. Studies of eddy dynamics 
and vertical fluxes also require measurements of local veloc- 
ities, along with passive and active scalars such as water 
vapor, temperature, and carbon dioxide. Although future 
remote sensors and certain current, but experimental, sys- 
tems may provide some of these measurements, it is clear 
that in situ observations remain necessary for most investi- 
gations. 

Combinations of remote and in situ measurements prom- 
ise to provide considerably more information than the sum of 
the parts. Many examples of combined observations are 
available in the literature, including Konrad and Robison 
[1973], Melfi et al. [1985], and Crum et al. , [1987]. 

This paper describes an experiment where the University 
of Wisconsin volume imaging lidar is used to portray the 
convective structure surrounding the flight path of the Ca- 
nadian National Aeronautical Establishment (NAE) Twin 
Otter aircraft. A unique feature of this study consists of 
precisely locating the aircraft measurements in the three- 
dimensional volume observed by the lidar. This allows point 
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TABLE 1. University of Wisconsin Volume Imaging Lidar 
Specifications 


Instructor 

Specification 

Value 

Transmitter 

laser 

Nd : YAG 


average power 

20 W 


repetition rate 

30 Hz 


wavelength 

106.4 nm 


pulse duration 

6 ns 

Receiver 

diameter 

0.5 m 


angular scanning rate 

25° s _l 


APD quantum efficiency 

35% 


range resolution 

7.5 m 


optical bandwidth 

1 nm 


average data rate 

~0.5 Gb h" 1 

Data Processing/Storage 

Function 

DEC VAX 11/750 

system control and data storage 

DEC LSI 1 1 

data transfer and formating 


CSPI array processor real-time data processing 

Stardent GS-1000 real-time graphics 

2.6-Gb write-once optical disk data storage 


by point comparisons between the aircraft measurements 
and the lidar-observed structure. 

Instrumentation 
The Volume-Imaging Lidar 

The volume-imaging lidar (VIL) is an elastic backscatter 
lidar designed to image the four-dimensional structure of the 
atmosphere. This system couples an energetic, high pulse 
repetition rate laser with a sensitive receiver and a fast 
computer controlled angular scanning system. High- 
bandwidth data acquisition is sustained during extended 
experiments by using a 2.6-Gb write-once optical disk for 
data storage. A Stardent GS-1000 graphics computer pro- 
vides 1280 x 1024 pixel resolution lidar images with 16-bit 
pseudocolor for real-time control of data acquisition and 
data analysis. System specifications and a block diagram are 
provided in Table 1 and Figure 1 . 

The high sensitivity of the VIL allows observations of 
inhomogeneities in natural aerosol backscatter that reveal 
clear air convective structure. Convective plumes carry 
surface aerosols aloft into cleaner air. Since many of the 
aerosol particles are hygroscopic, the particles also increase 
in size as the relative humidity increases. Adiabatic cooling 
of expanding air in a rising parcel increases the relative 
humidity and thus increases the aerosol scattering cross 
section. As a result the lidar easily delineates the boundary 
between mixed layer air and the free atmosphere above. 
Parcels of dryer, relatively clean air entrained into the 
boundary layer at the base of the inversion are also easily 
detected as they are mixed downward into the convective 
layer. 

Aircraft Instrumentation 

The Twin Otter contains a Rosemount 8585 five-hole 
probe to provide three orthogonal components of atmo- 
spheric motion over a frequency range of 0 to 5 Hz. Aircraft 


motions were determined with three separate systems: (1) 
three-axis Decca Doppler radar coupled with a NAE- 
assembled package of accelerometers and rate gyros, (2) 
Litton-90 inertial reference system, and (3) an ARNAV 
model R-40-AVA-100 Loran-C navigation system. Temper- 
ature was measured with a Rosemount fast response 
102DJ1CG heated probe mounted on the port side of the 
aircraft nose and by a second, identical probe mounted on 
the nose boom fairing. Fast-response humidity and carbon 
dioxide measurements used the ESRI infrared gas analyzer 
developed by Agriculture Canada. These measurements 
were made in a high flow rate (—300 L -1 ) duct passing 
through the aircraft cabin. Slow-response carbon dioxide 
measurements were made with a LI-COR LI-6251 carbon 
dioxide analyzer and slow-response water vapor measure- 
ments were obtained using an EG & G model 137 Cambridge 
dew point sensor. 

A more complete description of the aircraft instrumenta- 
tion is given by MacPherson [1990]. 

Data Acquisition 

The data described in this paper were acquired on August 
3, 1989, as part of the NASA First International Satellite 
Land Surface Climatology Project (ISLSCP) Field Experi- 
ment (FIFE). The lidar was operated south of Manhattan, 
Kansas, near the intersection of Route 177 and Interstate 
Highway 70 at FIFE site 123 (3639-L1W). 

For the observations reported in this paper, the lidar was 
programed to repeatedly scan an atmospheric volume con- 
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Fig. 1. Schematic diagram of the volume-imaging lidar. 
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sisting of the elevation angles between 0.3° and 10° and 
between azimuth angles of 215° and 255°. The scan pattern 
consisted of a sequence of 80 elevation scans separated by 
0.5° in azimuth; lidar profiles were acquired at elevation 
angle increments of approximately 0.25°. It took 187 s to 
scan the volume completely. Each scan included approxi- 
mately 5600 lidar profiles, each of which consisted of 1024 
data points obtained at 15-m intervals. 

The Twin Otter acquired data along constant pressure 
altitude flight legs, which began directly over the lidar and 
extended approximately 20 km from the lidar. Real-time 
lidar displays of the convective boundary layer were used to 
measure boundary layer depth as the aircraft approached the 
lidar. These measurements were used by the lidar operators 
to select the flight altitude of the aircraft for the data 
segment. The altitude was normally selected at a fixed 
fraction of the mixed layer depth. The flight heading was 
chosen along the center of the lidar scan, unless lidar 
imagery suggested a reason for an alternate heading. The 
flight altitude was converted to geometric altitude above 
mean sea level and the altitude and heading were radioed to 
the aircraft in time for the pilot to stabilize altitude, heading, 
and airspeed before beginning the data segment. This pro- 
cedure increased the efficiency of data collection by posi- 
tioning the aircraft at the proper altitude even in the presence 
of changes in boundary layer depth. During one observation 
period, convective “streeting” was observed in the lidar 
images. In this case, the lidar operators were able to select a 
flight heading perpendicular to the “streets” to clearly 
delineate the roll pattern in the aircraft data record. 

Data Analysis and Navigation 

Superposition of aircraft measurements on the lidar imag- 
ery requires careful attention to aircraft navigation and to the 
advection of atmospheric structure by the wind. The size of 
clear air convective plumes is comparable to the depth of the 
convective boundary layer; thus to maintain a 10% position- 
ing error, the aircraft location must be known to an accuracy 
of 100 m or better in a typical afternoon boundary layer with 
a 1-km depth. 

At the start of each flight leg the aircraft crew triggered an 
event flag in the data record at the instant they judged that 
the airplane was directly over the lidar. This event flag was 
used to check the accuracy of the aircraft Loran position by 
comparing the recorded positions with the known position of 
the lidar. Since the lowest flight was flown at an altitude of 
— 100 m and ground observations indicated that the pilot was 
successful in flying directly over the lidar, the event flags are 
believed to have been triggered when the aircraft was within 
100 m of the lidar position. On the lowest-level flight the 
Loran showed the aircraft to be 565 m east and 796 m north 
of the lidar when the event marker was triggered. Because 
the experiment was located far from existing Loran trans- 
mitter chains, inside the “midcontinent gap,” the rather 
large position offset is not surprising. This displacement was 
assumed to result from a systematic error in the Loran 
position; thus this displacement was therefore subtracted 
from all Loran positions to form a corrected position. To 
gain a rough check on this correction, it was compared with 
positions recorded at the time of the event flag on other 
overflights. Corrected positions at the time of the event flags 
for flights at altitudes of —400 and —800 m were 158 m east, 
92 m south and 187 m north, and 37 m east, respectively. 


The 400-m flight began 15 min before the low-level flight 
and the 800-m flight began 14 min later. The position 
deviations at the time of the event flags on the higher flights 
seem consistent with an expectation for errors —100 m, 
considering the difficulty of piloting the airplane over the 
lidar and determining when the plane was directly above the 
lidar at the higher flight altitudes. 

The lidar position was determined from a U.S. Geological 
Survey topographical map, using measured distances be- 
tween the lidar van and the buildings pictured on the map. 
The error in this position is expected to be less than 10 m: 
this error does not influence the intercomparisons because 
the aircraft positions were determined relative to the lidar. 
The alignment of the lidar’ s azimuth with geographical 
coordinates was initially performed using a magnetic com- 
pass. Lidar returns from the top of a radio tower at a range 
of 12.9 km were later used, along with a map, to provide an 
exact value for the conversion from lidar azimuths to geo- 
graphical directions. 

Knowing the aircraft position with respect to the lidar 
solves only a portion of the navigation problem. The wind 
moves the atmospheric structure during the time it takes the 
aircraft to fly through the lidar-observed volume and during 
the time it takes the lidar to complete a volume scan. At —45 
ms -1 it required —300 s for the Twin Otter to fly from the 
lidar van to the end of the 15-km range observed by the lidar. 
Similarly, the lidar required —200 s to scan a complete 
volume. Since the wind speed was —10 m s _1 when these 
data were acquired, the atmosphere moved —3.3 km during 
the flight leg and —2.2 km during the time it took for one lidar 
scan. At least once during each flight leg the lidar scanned 
directly past the aircraft; at this point and only at this point 
could the aircraft data and lidar image be compared directly 
without compensation for wind displacements. At all other 
points along the aircraft flight leg the air parcel sampled by 
the aircraft appears at some other location in the lidar image 
because of the wind. Both the lidar images and the apparent 
flight path of the aircraft were transformed to correct for 
wind displacements. 

Lidar Image Computations 

Individual lidar returns were first normalized for laser 
energy fluctuations between laser firings and corrected for 
the inverse range squared dependence contained in the lidar 
equation. A correction was also applied for optical attenua- 
tion, using the approximation that the attenuation cross 
section was independent of position in the volume. These 
corrections are more fully explained by Schols and Eloranta 
[this issue]. All images displayed in this paper have been 
corrected for the distortion produced by the mean wind 
during the time required for the lidar to scan the volume. In 
the discussion that follows, all times are measured from zero 
at the time the aircraft flew over the lidar. 

Constant altitude plan position indicator (CAPPI) images 
were generated as follows: 

1 . A synthetic lidar profile was generated for each of the 
80 range height indicator (RHI) scans making up a volume 
image. This profile reproduced the corrected lidar return, 
which would be observed by a horizontal pointing lidar fired 
at the azimuth of the individual RHI scan if it were located at 
the altitude of the CAPPI display. At each range, data points 
in the horizontal pointing synthetic profile were created from 
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Fig. 2. Preparation of constant altitude plan position indicator 
(CAPPI) lidar images from volume scans. Profiles corresponding to 
a horizontal pointing lidar at the altitude of the CAPPI are generated 
from each range height indicator (RHI) scan. The data points 
between ranges 1 and 2 in the horizontal profile at altitude, h , are 
generated by linear interpolation between successive lidar profiles 
obtained along directions p and q . This process is repeated for all 
other pairs of adjacent profiles in the RHI scan to complete the 
horizontal profile. Horizontal profiles from all RHI scans in the 
volume scan are then used to compute the CAPPI image. 


the corrected data points in the original RHI scan, as shown 
in Figure 2. The synthetic data points are computed from a 
linear interpolation between the data points in the nearest 
two lidar profiles. 

2. Assuming a constant wind velocity, the individual 
profiles are translated to the position the air occupied at the 
instant when the aircraft passed over the lidar van. This 
process corrects for the distortion introduced by the wind in 
the -200 s required to complete a volume scan. Notice that 
since the profiles that make up the CAPPI are synthesized 
from many individual profiles in the original RHI scan, this 
procedure assumes that the CAPPI profile was acquired 
instantaneously and neglects the 1- to 2-s time delay between 
the first and last shot in the RHI. 

3. Finally, the missing pixels between the synthetic 
profiles are filled in by linear interpolation. This produces a 
complete CAPPI image depicting the atmospheric structure 
at the CAPPI altitude at the instant the aircraft passed over 
the lidar. 

To superimpose the aircraft-observed data points on the 
wind-corrected CAPPI image, the position of each aircraft 
data point must be translated back to the position occupied 
by the air parcel at the instant when the aircraft passed over 
the lidar. The vector displacement for each point is the 
negative of the mean wind vector, multiplied by the time 
between the measurement and the time the aircraft passed 
over the lidar. With this displacement the aircraft data can be 
plotted directly on the CAPPI image. This procedure cor- 
rects for the geometric distortion introduced by the mean 
wind field; however, it does not account for temporal evo- 
lution in the structure. The comparison between aircraft 
observations and lidar imagery can be expected to degrade 
as the time separation between the lidar and the aircraft 
observations of a particular air parcel increases. Superposi- 
tion of aircraft data on lidar RHI images is slightly more 
difficult. Individual RHI images are acquired rapidly (1-2 s) 
so that there is little wind induced distortion in an image; 
however, because the aircraft requires —300 s to complete 
the flight leg, a composite RHI image must be formed that 


portrays the air parcels actually penetrated by the aircraft. In 
the following discussion, times are measured from the time 
at which the aircraft passed directly over the lidar van. 

The composite RHI images were formed as follows: 

1. The flight track of the aircraft over the ground was 
plotted and determined to be very close to a straight line. 

2. A straight line was fitted to the ground path of the 
aircraft and the average speed of the aircraft along the path 
was calculated. 

3. Each aircraft position on the best fit line was displaced 
upwind by the vector wind multiplied by the time. This new 
wind corrected straight line of positions indicated where the 
aircraft-sampled air parcels were at the time where the 
aircraft passed over the lidar. 

4. For each RHI scan a pair of radial lines 0.25° in 
azimuth above and below the RHI scan were constructed. 
These lines were also displaced upwind a vector distance 
equal to the vector wind multiplied by time. 

5. The points at which these lines intersected with the 
wind-corrected aircraft position line were then computed. 
These intersection points indicated the upper and lower 
range limit in the RHI image where this RHI azimuth 
provided the best match in space and time to the air parcel 
sampled by the aircraft. 

6. Image points from the RHI scan between these range 
limits were then mapped onto the composite image. 

7. All RHI scans from this volume scan were checked for 
intersecting segments and a composite image was formed. 

8. The aircraft data were superimposed onto the com- 
posite RHI image, with each data point plotted on the 
straight-line fit to the aircraft path at the position nearest to 
the aircraft position at the time the data point was acquired. 

Three-dimensional contour surfaces of lidar backscatter 
were generated as follows: 

1. The position of each profile in the lidar scan was 
displaced upwind by a distance equal to the wind speed 
multiplied by the time. 

2. Inside the transformed data volume a contour algo- 
rithm computed elemental polygons separating data points 
whose lidar backscatter exceeded a selected contour value 
from all other points. 

3. Perspective transforms and a Gouraud shading algo- 
rithm were used to display the contour surface as a solid 
body image. 

4. A system probatoire pour 1’ observation de la terre 
(SPOT) panchromatic satellite image of the terrain under the 
lidar-scanned volume was displayed in perspective to show 
the local topography. 

5. Each point in the aircraft flight path was displaced 
upwind by a distance equal to the wind speed multiplied by 
the time. This corrected flight path was then displayed as a 
line penetrating the lidar image. 

Data Presentation 

This paper presents data measured on August 3, 1989. The 
sky was partly covered with fair weather cumulus clouds; 
cloud coverage is shown in Plate 1. At the time of the 
measurement the mean boundary layer depth was —850 m, 
with the tops of the highest cumulus observed at —1300 m. 
Cloud base altitude was —910 m. All heights are expressed 
as the distances above the lidar. 

Wind profiles measured with the lidar and averaged over 
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Plate 1. Three-dimensional solid contour image of clouds and the aircraft flight path derived from the lidar scan 
obtained between 1855:43 and 1858:36 UT. This image has been corrected to represent the position of the clouds when 
the aircraft passed over the lidar. Wind corrections used the lidar-measured wind of speed of 10.2 m s _l and direction 
of 215.2°. Since the lidar penetrates only a short distance into these clouds, the image shows the cloud base and 
unobscured portions of the cloud facing the lidar. The terrain image was computed from a panchromatic SPOT satellite 
image provided by the SPOT Image Corporation. (Copyright CNES 1987.) 


the —70 km squared area of the scan were nearly indepen- 
dent of height between 200 m and 900 m, with a speed of — 10 
m s" 1 and a direction of —213°. Lidar wind profiles are 
presented in Figures 3 and 4. This data set was chosen for 
analysis because the strong winds provide a stringent test of 
wind corrections and navigation. However, because of me- 
chanical turbulence that mixes the aerosol structure, these 
rather windy conditions are not ideal for acquiring lidar 
images with good definition in the lower part of the boundary 
layer. 

Figure 5 presents variables measured by the Twin Otter on 
a flight leg flown at an altitude of 800 m above the lidar and 
beginning at 1854:53 UT (1354:53 CDT) and ending at 
1901:33 UT. Instantaneous fluxes are also plotted in Figure 
5. Primed quantities have been computed by subtracting the 
mean value over the flight leg; no other detrending or 
filtering has been applied. Distances along the flight path 
were taken from the Loran position, using the Loran offset 
measured on the lowest-altitude flight. Aircraft ground speed 
calculated from the Loran position was 46.0 ms -1 averaged 
over the duration of the data record. Aircraft altitude was 


computed from a combination of radar altimeter and pres- 
sure altitude data. The radar altimeter has a maximum range 
of 779 m and thus could not be used to directly measure the 
altitude of this flight leg. The radio altimeter was used to 
measure altitudes of 436, 436, and 102 m above the lidar on 
three flight legs which began at 1817, 1826, and 1841 UT, 
respectively. The difference in pressure altitudes between 
these overpasses and the pressure altitude of the flight 
beginning at 1855 was corrected for the difference in temper- 
ature of the standard atmosphere and virtual temperature 
measured by the aircraft to generate true altitude. These 
three computations produced nearly identical results of 803, 
802, and 803 m above the lidar for the 1855 leg. 

Wind measurements were acquired both from aircraft 
instruments and from lidar observations. Lidar winds were 
computed by observing the drift of naturally occurring 
inhomogeneities in the atmospheric aerosol content using 
the algorithm described by Schols and Eloranta [this issue]. 
The measurements represent an average over the area of the 
lidar scan (—70 km 2 ) and the time interval between the scans 
(187 s). Two lidar wind profiles were measured during the 
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Fig. 3. Lidar-measured wind speed profiles for the periods 
1855:43 to 1858:36 UT (solid squares) and 1858:50 to 1901:42 UT 
(open squares). Leg-averaged wind speeds computed from the 
aircraft data are also presented. Aircraft wind computations are 
shown using both the Litton-90 INS (solid circle) and the Loran 
(open circle) to determine aircraft motion. 


aircraft flight leg (Figures 3 and 4). Aircraft winds represent 
the wind averaged over a 15-km data segment flown through 
the lidar scan. Computation of winds from the aircraft 
observations required knowledge of the aircraft velocity 
averaged over the time of the flight leg. This information was 



Fig. 4. Lidar-measured wind direction profiles for the periods 
1855:43 to 1858:36 UT (solid squares) and 1858:50 to 1901:42 UT 
(open squares). Leg-averaged wind speeds computed from the 
aircraft data are also presented. Aircraft wind computations are 
shown using both the Litton-90 INS (solid circle) and the Loran 
(open circle) to determine aircraft motion. 



Fig. 5. Airborne measurements showing fluctuations about the 
flight leg mean value for vertical wind, temperature, humidity, and 
carbon dioxide concentration fluctuations as a function of range 
from the lidar. Data were acquired at a mean altitude of 803 m above 
the lidar in the time interval between 1854:53 and 1901:33 UT. No 
filtering or detrending has been applied to these records. Flux 
quantities are simple products of the displayed vertical velocity and 
scalar fluctuation quantities. 


derived from the Loran positions at the ends of the data 
segment and the duration of the flight. A comparison of the 
Loran positions with the Litton inertial navigation system 
(INS) showed a 0.485 m s _l drift of the INS at 304° with 
respect to the Loran positions over the duration of the data 
segment. 

The wind measured by the lidar at an altitude of 800 m was 
10.2 m s** 1 at 216.3°; this was obtained by averaging the two 
lidar measurements obtained during the aircraft data run. A 
sequence of seven lidar wind measurements obtained be- 
tween 1824 UT and 1901 UT, at an altitude of 800 m, shows 
speeds between 10.1 m s -1 and 10.9 m s _I , with directions 
between 213° and 217°. Aircraft wind measurements show 
speeds of 10.32 m s _l at 216.2°, using the Loran, and 10.35 
m s" 1 at 218.9°, when using the Litton-90 INS to determine 
aircraft velocity. A comparison of the lidar winds computed 
at the time of the aircraft flight leg and the aircraft winds are 
provided in Figures 3 and 4. Note the excellent agreement 
between measurements; at the airplane flight altitude the 
successive lidar measured speeds agree to within 0.1ms" 1 
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Plate 2. A composite RHI scan where aerosol structure in the volume scan is displaced to correct for motion caused 
by the mean wind between the time the aircraft sampled an air parcel and the time the lidar scanned the same parcel. 
The time when the volume scan was completed is displayed on the composite image. Water vapor mixing ratio 
variations for the flight beginning at 1854:53 UT on August 3, 1989, are superimposed on the lidar image. The distance 
of the aircraft from the lidar is shown on horizontal axis. An airplane symbol marks the location where aircraft and lidar 
measurements are coincident. 


and direction agrees to within 2°. All lidar and aircraft speed 
measurements agree within 0.2 m s~ ] and directions agree 
within 4°. 

The ability of navigation and wind corrections to properly 
place aircraft measurements on the lidar images of air 
parcels can be judged by superimposing aircraft water vapor 
measurements on lidar images. Water vapor is selected 
because relative humidity strongly influences the optical 
scattering cross section of naturally occurring aerosols. 
Many of the aerosols are composed of hygroscopic materials 
[Fitzgerald and Hoppel , 1984; Sale mink et al 1984], which 
grow by the accretion of water when the relative humidity 
exceeds the deliquescent point of the aerosol. This flight leg 
was just beneath the bases of small fair weather cumulus 
clouds, the relative humidity was large in the convective 
plumes and relatively low in the parcels entrained from 
above the mixed layer. These conditions produce a sharply 
defined lidar image to compare with the humidity measure- 
ments. 

A composite RHI lidar scan, with a plot of the water vapor 


mixing ratio superimposed, is displayed in Plate 2. Wind 
corrections to the lidar image and the aircraft flight path were 
derived using the winds measured with the lidar (10.2 m s _l 
at 215.2°). 

Aerosol structure in the composite RHI is highly corre- 
lated with the aircraft-based water vapor measurements: 
regions of enhanced backscatter correspond to high values 
of humidity. In particular, note the tongue of dry air at a 
range of —7 km. A close correspondence between the 
brightness of small-scale aerosol structure and the humidity 
trace is seen at all ranges between 5 and 9 km. The lidar 
scanned near the aircraft at a range of —6.8 km, providing 
coincident measurements at this range. Thus temporal evo- 
lution of the structure does little to degrade the aircraft lidar 
comparison at the left-hand side of Plate 2. Plate 3 shows the 
RHI composite generated from the lidar volume scan ob- 
tained 187 s after Plate 2. In Plate 3 the aircraft and lidar data 
are nearly coincident at the extreme range of the lidar image. 
When Plate 2 is compared to Plate 3, temporal evolution of 
the structure is clearly evident. In Plate 3 the comparison 
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Plate 3. The composite RHI scan computed from a lidar volume scan which began 187 s after the scan used to produce 
Plate 2. In other respects this figure is identical to Plate 2. 


between lidar and aircraft observations are best beyond 10 CAPPI observed 187 s later (Plate 5). In this image the large 

km. Note that the tongue of dry air shown near 7 km in Plate dry incursion of air that appears under the aircraft path 5 km 

2 has disappeared in Plate 3. The entire region between 10 south of the lidar is partly filled in by new plumes: one of 
km and 12.5 km has also become darker in Plate 3. This is these was just visible on the previous image (4.5 km south, 
consistent with the downdraft region apparent in the vertical 6.3 km west). In Plate 5 the image correlates well with the 
velocity record (Figure 5). humidity trace at the southwest end of the picture but not as 

The close correspondence between aerosol structure and well in the northeast. For example, note the dry region just 

humidity is further illustrated in Plate 4 where humidity is beyond 8 km south of the lidar: this now corresponds nicely 

plotted on a wind-corrected CAPPI. In the region of the to a dark spot in the image. 

picture near the aircraft the comparison is good: for exam- Figure 5 illustrates the well-known difficulty of measuring 
pie, the small plume under the aircraft path 4.5 km south of fluxes near the top of the mixed layer [Lenschow and 
the lidar appears as a sharp peak in the humidity trace. Up to Stankov , 1986]. In the course of the 15-km flight, only about 
a point 6 km south of the lidar, even very small features in five significant “flux events” are evident. Clearly, the sta- 

the image match with the details of the humidity trace, tistical reliability of an estimate made from such a low 

Beyond 6 km south, small discrepancies appear. The dry number of events is poor. An examination of the vertical 

region just beyond 6 km appears to correspond to a dark spot velocity trace shows the largest vertical velocity events near 

just north of the aircraft track, and a dry region just beyond the lidar. Fluxes computed using this record will be very 

8 km south appears to correspond to a dark spot south of the sensitive to the detail of any detrending or low-pass filtering 

aircraft track. This is caused by the time difference between of the vertical velocity trace. 

the aircraft observations and the lidar image. Temporal The “top hat” humidity plume measured by the aircraft 10 
evolution has modified the structure so as to reduce the km from the aircraft (Figure 5) marks the plume directly 

correlation. The extent of the evolution in structure can be under a fair weather cumulus cloud shown at ~7 km south, 

judged by comparing the CAPPI displayed in Plate 4 with the 10 km west in Plate 1. Reference to the vertical velocity 
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Plate 4. Wind-corrected CAPPI at the altitude of the aircraft flight measurements derived from the lidar volume 
scan obtained between 1855:43 UT and 1858:36 UT. Also plotted along the aircraft flight path is the variation of water 
vapor mixing ratio from its flight leg mean value. The image represents the position of aerosol features when the airplane 
passed over the lidar. Wind corrections were made using a mean wind speed of 10.2 ms -1 and direction of 215.2°. 


already penetrated significantly into the mixed layer. The 
next CAPPI (Plate 5) shows elements of mixed layer air 
rising into much of the region occupied by the downdraft in 
first image; evidently, this parcel was nearing the end of its 
downward motion when sampled by the aircraft. Thus even 
though this downdraft represents the largest single flux event 
observed along the flight leg, it is likely that the local fluxes 
in the parcel were significantly larger prior to the aircraft 
observation. At an earlier time, before the negative buoy- 
ancy of the plume had time to slow the downward motion, 
we would expect stronger downdraft velocities. 

Three of the prominent features in this record appear to 
have been sampled at times other than the time of maximum 
flux contribution. This indication that strong flux contribu- 
tions by an individual thermal are transient adds to the 
difficulty of obtaining reliable flux measurements in the 
upper part of the mixed layer. Regions of strong vertical 
motion are sampled infrequently. Regions contributing large 
fluxes as a result of a combination of both large vertical 
velocities and large excursions in scaler quantities are there- 
fore even more infrequent. 


inside the top hat shown in Figure 5 indicates that this plume 
is no longer supported by rising air. A second cloud at 9 km 
south and 13 km west in Plate 1 corresponds to the humidity 
top hat between 12.5 and 13.5 km in Figure 5. This cloud also 
shows little support from vertical velocity under the cloud. 
According to the classification of Stull [1985] both of these 
clouds are passive. It appears likely that if the aircraft had 
sampled fluxes under these clouds during their active phase, 
much larger fluxes would have been measured. Although the 
cloud shown between 12.5 and 13.5 km in Figure 5 exhibits 
little vertical velocity under the cloud, a flux contribution is 
provided by the downdraft directly to leeward of the cloud. 
The downdraft is clearly evident as a small dark spot in Plate 
4. 

Examination of the instantaneous fluxes of water, temper- 
ature, and carbon dioxide show the strongest single contri- 
bution from the downdraft located at —7 km from the lidar 
(see Figure 5). This downdraft exhibits large fluctuations in 
the scaler quantities but only rather a small downward 
velocity (—0.25 ms -1 ). It seems reasonable that the vertical 
velocity is small, since the negatively buoyant plume has 
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Plate 5. Wind-corrected CAPPI derived from the lidar volume scan obtained between 1858:50 and 1901:42 UT. This 
image is derived from the volume scan data acquired 187 s after Plate 4; otherwise, this figure is identical to Plate 4. 


Conclusions 

Volume-imaging lidar observations can be used to portray 
the atmospheric structure around the flight path of an 
instrumented aircraft. Because the current lidar technology 
requires approximately 3 min to map the three-dimensional 
structure, images must be corrected for distortion created by 
the wind. This paper shows that aircraft and lidar measure- 
ments can be successfully combined, even in the presence of 
10ms" 1 winds. Loran provided aircraft navigation informa- 
tion; however, it was necessary to subtract a constant 
displacement vector from the Loran positions. Because a 
low-level flight over the lidar van was available to provide 
this offset, aircraft data matched lidar imagery without any 
additional navigation corrections. In future experiments, 
aircraft location information is expected to improve greatly 
with the introduction of global positioning system naviga- 
tion. Evolution of the convective structures can be judged 
from successive lidar images; however, the rate of turbulent 
change is rapid, and higher lidar scan rates would be 
desirable. 

A sample flight leg illustrates the well-known difficulties of 
measuring fluxes in the upper part of the mixed layer. It is 


clear that flight legs much longer than the 15-km legs used in 
this experiment are required to provide stable estimates. 

Winds measured by lidar and aircraft show excellent 
agreement. When Loran was used to determine average 
aircraft velocity, flight-leg-averaged horizontal winds mea- 
sured by the aircraft and area-averaged winds measured by 
lidar agree to within 0.2 ms" 1 in speed and 1° in direction. 
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Abstract 


This paper presents automatic methods for obtaining convective boundary layer mean 
depths, cloud base altitudes, cloud top altitudes, cloud coverages, and cloud shadows 
using the University of Wisconsin Volume Imaging Lidar (VIL). These methods provide 
mesoscale observations representing the 70 km 2 scanning area of the VIL. All data obtained 
with the VIL during the 1989 FIFE field experiment are analyzed and the results axe 
summarized. The reliability of these methods is verified by comparing the results with 
radiosonde profiles, satellite images, and visual estimates from Range Height Indicator 
scans of the VIL. 
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1. Introduction 


Convective boundary layer (CBL) mean depths, cloud coverages, cloud base alti- 
tudes, and cloud top altitudes are traditionally measured with satellites, radiosondes, 
and weather stations. Most applications try to describe the atmospheric conditions in 
large areas prompting a need for area averaged measurements. Representative area av- 
eraged observations are generally achieved by applying time-averaging to ground-based 
point measurements. However, a long averaging time interval does not necessarily guar- 
antee a representative area average, since the point measurements can be biased by local 
variations. The accuracy of short averages is limited by convection scale inhomogeneities. 
This paper introduces methods for obtaining area averaged measurements of CBL mean 
depths, cloud coverages, cloud base altitudes, and cloud top altitudes with the University 
of Wisconsin Volume Imaging Lidar (VIL). These measurements provide stable averages 
over ~70 km 2 area with a 3- minute time resolution. 

The VIL [Eioranta and Forrest , 1992] is a scanning lidar capable of mapping elastic 
backscatter from aerosols in a large volume. The current VIL system configuration is 
summarized in Table 1. The VIL has a 10 MHz analog-to-digital conversion rate and 1 
kilobyte data buffer length that allow a 15 km detection range with a 15 m range resolution. 
The data acquisition system normalizes the backscatter profiles for laser output energy 
fluctuations, corrects the r tinge square dependence, and forms a natural logarithm of the 
signal before saving it in compressed form on optical disks. Scanning is performed with a 
beam steering unit mounted on the optical axis of the system. Figure 1 illustrates a VIL 
volume scan. For the convective boundary layer studies, the VIL typically scans 15° high 
and 40° wide sectors with 0.33° and 0.5° angular step sizes, respectively. The VIL devotes 
about 3 minutes to scan this ~100 km 3 volume. 
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The VIL was operated in the First ISLSCP (International Satellite Land Surface 
Climatology Project) Field Experiment (FIFE) site [Se/lers et a/., 1992] from July 26 to 
August 11, 1989. It recorded ~ 2000 boundary layer volume scans consisting of ~ 20 
gigabytes of raw data. In earlier studies [Boers and Eloranta, 1986; Crum et ah, 1987; 
Nelson et a 1., 1989], the convective boundary layer mean depth and entrainment zone 
thickness were measured by a visual inspection of Range Height Indicator (RHI) scans. Due 
to the enormous amount of raw data from the 1989 FIFE program, the visual inspection 
was impractical. Therefore, automatic analysis methods for measuring CBL mean depths, 
cloud base altitudes, cloud top altitudes, cloud coverages, and cloud shadows were essential 
in order to analyze all data from the experiment. 
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2. Methods 


2.1. Convective Boundary Layer Mean Depth in Clear Air 


The fair-weather convective boundary layer is characterized by rising thermals of 
aerosol rich air and sinking streams of air from the free atmosphere. These vertical motions 
cause irregularities in the boundary layer top which can be observed from the RHI scans 
of the VIL (see Figure 2). Because of insufficient sampling of the undulating CBL top, 
traditional point measurements of the CBL depth axe not representative. 

The CBL top is typically marked with a sudden lidar backscatter intensity change 
due to different aerosol concentrations in and above the boundary layer. When visually 
inspecting RHI scans, we have adopted the definition by Deardorff et al. (1980) for the 
CBL mean depth: the altitude which has equal areas of free air below and boundary layer 
aerosols above. In automatic estimation, we determine the CBL mean depth from the 
horizontal signal variance [Hooper and Eioranta, 1986]. 

The horizontal aerosol variations are intense at the altitude of the CBL mean depth 
because of undulations at the top of the CBL. To calculate the vertical profile of horizontal 
signal variations, each shot in a volume scan is first high pass filtered with an 1 km 
long median filter. This eliminates atmospheric extinction and scattering from uniformly 
distributed aerosols. Using a 20 m vertical spacing, backscatter signals from each altitude 
range are mapped on horizontal planes. Then, the variance of the aerosol signal on each 
horizontal plane is calculated. The lowest altitude local maximum peak of the variance 
profile marks the convective boundary layer mean depth. 
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The bottom-up approach eliminates errors due to strong aerosol variability caused 
by aerosol layers above the CBL. Local maximas induced by dust clouds are rejected by 
requiring that the peak is stronger than the average of the signal variance profile. Random 
variance fluctuations due to signal noise can disturb detection of maximum variance when 
the difference in backscatter from boundary layer aerosols and free air is small. To prevent 
selection of a random local maximum we check not only that the points either sides of 
the maximum point are smaller than the maximum, but also that these adjacent points 
are larger than the next closest points, i.e. <r(z max — 2 A z) < cr(z max — A z) < cr(z max ) > 
v(zmax + A z) > cr(z max + 2A z). 

The automatic method provides a convective boundary layer mean depth estimate 
representing an average over ~ 70 km 2 . Figure 2 shows an RHI scan, the VIL signal 
variance profile, and the estimated CBL mean depth. The automatic method is tested 
by comparing the results with manually inspected estimates. The manual estimation is 
performed by determining the altitude which has equal areas of boundary layer aerosols 
above and sinking air below. Figure 3 compares manually and automatically determined 
convective boundary layer mean depths and estimates from radiosonde-based potential 
temperature profiles on July 28, 1989, from 8:00 to 10:00 CDT. Because of the large 
number of RHI scans during the session, the manual inspection is performed for three RHI 
scans in the same volume scan every fifth volume scan. The variability of manual estimates 
increase in time due to growing undulations at the top of the CBL. The automatic estimates 
are consistent and show a smooth CBL growth over time. The fluctuation in the automatic 
results at 9:50 CDT is due to a missing data block in the volume scan; the transmitter 
laser beam was turned off when an aircraft flew into the field of view. 

The plot of potential temperature profiles used for estimating the CBL mean depths 
in the previous case is plotted in Figure 4. Traditionally, the CBL top is located from a 
jump in the potential temperature profile. Since the radiosonde is a point measurement 
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system, the potential temperature soundings are sensitive to the spatial distribution of the 
thermal plumes. The potential temperature profile can be measured through a rising plume 
or through a downward stream of free air penetrating deep in the CBL. This can lead to 
uncertainties which are in the order of the CBL mean depth. The potential temperature 
profile at 8:01 CDT shows a 100-150 m CBL mean depth, which is consistent with the 
VIL measurements. The CBL depth is very difficult to determine from the profile at 9:36 
CDT. The first potential temperature jump is at ~ 300 m, which is ~ 200 m lower than 
the VIL estimates. The profile also shows an instability at ~ 700 m, which was probably 
due to the radiosonde rising through a thermal. Thus, the CBL top is most probably at 
~ 730 m, where the potential temperature starts to rise again. However, the radiosonde 
profile was most probably measured through a high plume top. 

2.2. Convective Boundary Layer Mean Depth with Clouds 

When cloudiness increases, echoes from clouds can bias the CBL depth estimate by 
dominating the variance. In such cases, the maximum variance represents the altitude of 
the maximum cloud echo variability instead of the convective boundary layer mean depth. 
The method for determining CBL mean depths in the presence of a few low altitude clouds 
is similar to the clear air case described in the previous section, except the shadowed regions 
behind clouds are removed before computing the horizontal variance. 

Figure 5 presents a comparison between manually and automatically determined CBL 
mean depths from the 1989 FIFE program. The manual and automatic estimates correlate 
well. Only data with the cloud coverage less than 10% is plotted, since in presence of more 
clouds, both automatic and manual method can give biased results due to cloud shadows. 
When the CBL mean depth is below 200 m, both methods may give statistically poor 
estimates, since the lidar images often shows only the tops of the plumes. Therefore, 
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there may not be enough information to determine the CBL mean depth. Otherwise, the 
deviations are mostly due to under-sampling of undulating CBL top in the manual method, 
which causes typically 100-200 m variations in estimates. In order to understand better 
the scatter in this plot, some of the most extreme discrepancies are reanalyzed. Boxes a 
and c mark CBL top estimates which were difficult to determine with the manual method 
due to low contrast between boundary layer aerosols and free air. Box b marks points, 
where a haze layer just above the CBL was manually misinterpreted as the CBL top. After 
reinterpretation, the manual estimates in marked boxes are within 100 m of the automatic 
measurements. 

In the presence of a significant convective clouds coverage, the definition of a boundary 
layer mean depth is unclear and the previously presented convective boundary layer mean 
depth determination is invalid. For example, near a thunderstorm cloud with a cloud top 
at 10 km, the local boundary layer would be 5 km by the definition we have adopted. 

In the presence of low clouds, the cloud base and cloud top altitudes are presented. In 
these cases, the values of the mean depth must be interpreted with a great caution. Figure 
6 presents an RHI scan along with a signal variance profile and the CBL mean depth 
estimate. Cloud echoes cause a strong variance. The automatically determined CBL mean 
depth is slightly higher than the manual estimate. 

I 

I 

2.3. Cloud Base Altitude 

The lidar backscatter echo from a cloud is detected as a strong positive signal gradient 
(see Figure 7). The first signal gradient equal to the maximum slew rate of the detector 
amplifier is used to define the cloud base. The gradient is calculated from three subsequent 
data points to prevent misinterpreting noise induced signal gradients as clouds. To obtain 
representative cloud base estimates, all shots in a volume scan axe inspected for cloud 


7 



echoes and the lowest cloud echo defines the cloud base altitude. This method provides 
the lowest cloud base altitude in the ~ 70 km 2 VIL scanning area. Sometimes the VIL 
backscatter signal includes echoes from ground objects and dust clouds, which are similar to 
cloud echoes. Therefore, every cloud echo below 200 m is omitted in the cloud calculations. 

Figure 8 presents an RHI scan with a histogram of the first cloud return altitudes and 
an estimated cloud base. The cloud base is measured with ~50 m accuracy. Clouds are 
good scatterers and are clearly visible in an RHI scan. However, when inspecting cloud 
bases from the RHI scans, the visual difference between thermals and clouds depends 
on the color enhancement selected for the display terminal. For example, some clouds in 
Figure 8 seem to extend below the detected cloud base, since the chosen color enhancement 
shows backscattering from aerosols and clouds with similar grey shades. An accurate visual 
estimation of the cloud base from the RHI scans requires a well trained eye and careful 
setting of enhancement levels. 

The automatic cloud base altitude estimates correlate well with the manually deter- 
mined estimates. Figure 9 represents the manual and automatic cloud base estimates for 
the 1989 FIFE experiment data. Only those cases with cloud coverage greater than 5% 
axe analyzed. To reduce statistical fluctuations from sampling differences, a 15-minute 
minimum filter is applied to the results. Most manual estimates are greater than or equal 
to automatic estimates. The manual method samples only ~1 % of the RHI scans and it 
is unlikely to sample lowest clouds in the volume scan causing systematic over estimation. 
We expect to have statistical fluctuations of ~50 m, since the vertical resolution is limited 
by the elevation angle step size used in the volume scan. Boxes a— d mark regions where 
manual estimates are substantially higher than the automatic measurements. The visual 
estimates for these cases have been re-examined to explain the differences. One sample in 
box c is due to a dust cloud rising up to 500 m, which caused the automatic method to fail. 
The other measurements in box c were made when there was a cloud layer above a few 
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fair-weather cumulus clouds. In this case, the manual method missed the lower clouds due 
to insufficient sampling. Boxes e and f mark manual estimates which are more than 200 m 
lower than the^automatic measurements. These are due to choosing an inappropriate color 
enhancement in visual inspection of the RHI scans showing thermal plumes and aerosol 
swelling as clouds causing under estimated manual cloud base altitudes. After reinterpre- 
tation of the RHI scans, the manual estimates in the marked boxes return within 50 m of 
the automatic estimates. When values in boxes a— f are reanalyzed, the root-mean-square 
difference between the the automatic and manual results is ~67 m. 


2.4. Cloud Top Altitude 


Clear lines of sight are required for a reliable cloud top altitude estimate, since the 
extinction in a cloud often blocks the lidar signal before the cloud top is reached. A lidar 
pointing directly to zenith can only detect tops of optically thin clouds. However, the 
three-dimensional scanning capability of the VIL extends the visibility of cloud tops by 
looking through holes in the cloud deck. This makes the cloud top detection possible for 
scattered, optically thick clouds. 

Assuming that 1) the horizontal scales of clouds are approximately equal to CBL 
mean depth, 2) the shape of clouds is square, and 3) the clouds are uniformly distributed, 
we can derive am approximate equation for the maximum detectable cloud height h as a 
function of the observation angle, fractional cloud coverage C, CBL mean depth z,-, and 
penetration depth d p of the lidar signed into a cloud as 


h 


{{is- x ) Zi+i r 


tan# 


( 1 ) 


where 6 is the maximum observation angle (< 90°). The penetration depth of the VIL 
signed into clouds is generedly greater than 200 m. This equation provides rough limits of 
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the cloud top detection in presence of optically thick clouds. For example, if z, = 1 km, 
C = 0.25, and 6 = 10°, then the cloud tops below ~1200 m can be detected, even for 
optically thick clouds. 

The maximum altitude of the last returns from clouds in a volume scan defines the 
cloud top altitude. The last cloud return is marked by a rapidly falling signal after cloud 
maximum return (see Figure 7). Because a strong extinction in an optically thick cloud 
may cause a similar signature, the rest of the backscatter profile is checked for remaining 
backscatter. The ratio of lidax profiles penetrating through the cloud layer to all cloud 
top returns indicates the reliability of the cloud top estimate. The cloud top altitude is 
manually detected from an RHI scan as the highest visible cloud structure. Figure 10 
shows an RHI with a histogram of the last returns from clouds and a cloud top estimate. 
The cloud top estimate is higher than one would estimate from this RHI scan, since the 
automatic cloud top is calculated using the whole volume scan. 

The automatic cloud top altitude estimates correlate with the manual estimates, as 
shown in Figure 11. A 15-minute maximum filter is applied to the results to reduce 
statistical fluctuations due to sampling rate differences between the methods. Box a marks 
the point where a dust cloud rising up to 500 xn caused the automatic method to fail. Boxes 
b and e mark points where a human error caused spurious values in the manual estimates. 
Boxes c, d, and f mark points where a higher cloud layer was missed by the manual 
method, since the vertical range of the display was set too small. After reinterpretation 
of RHI scans, the manual results marked in boxes b— f are within 50 m of the automatic 
estimates. In cases of partly transparent clouds, the manual method sometimes provides 
higher cloud top estimates, since it can detect faint cloud echoes above the top of the 
first cloud layer. When clouds are optically thick, both methods have the same visibility 
problems. The automatic method, however, samples each shot and thus provides more 
representative estimates due to better sampling coverage. 
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2.5. Cloud Coverage and Cloud Shadows 


Cloud coverage is defined as the percentage of ground surface area with clouds directly 
above. The cloud echoes and signal echoes are mapped into a horizontal grid according to 
their horizontal positions. The ratio of grid points with cloud echoes to all grid points with 
signal echoes defines the cloud coverage percentage. The VIL can map clouds on areas as 
small as ~ 2000 m 2 within the scanning area. Cloud coverage estimation is also possible 
using a larger portion of the scanning area up to ~ 70 km 2 . 

The measurement geometry limits the accuracy of the cloud coverage estimates. The 
coverage decreases towards the end of the scan range since the lowest cloud elements at 
close range may block the view of clouds with slightly higher bases. To minimize this 
error, the cloud coverage is calculated using 5 km long and 5 km wide grids starting from 
the closest possible range for the cloud base altitude and the maximum elevation angle of 
the scan. In terms of elevation angles, this means that the coverage of a cloud layer at 1 
km is measured using scanning angles between 6.5° and 15° to enable complete cloud base 
coverage when the lowest cloud elements are less than 100 m below the highest cloud base. 
Figure 12 shows the cloud echo projections to the ground and the spatial coverage of the 
backscattering returns. Cloud coverage is about 25% in this case. 

To verify the cloud mapping capability of the VIL, cloud coverage estimates axe com- 
pared with satellite imagery data. Figure 13 compares cloud mapping based on VIL data 
with the visual wavelength band (0.572-0.698 fim) image of the satellite-based Advanced 
Very High Resolution Radiometer (AVHRR) instrument on the NO A A- 11 satellite plat- 
form. The cloud base altitude was 1.4 km; the cloud top altitude was 2.4 km. The view 
angle of the AVHRR was 21.9° from zenith and 258.5° clockwise from North. Similar 
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cloud structures can be seen in both images. The main differences are due to the spatial 
resolution differences; the areal resolution of VIL is about five hundred times better than 
the 1 km pixel resolution of the AVHRR. 

The percentage of ground shadowed by clouds generally differs from the cloud cov- 
erage. In most studies requiring cloud shadow information, shadows are estimated using 
cloud coverage measurements with assumed cloud shape. This works well in large areas, 
but when large scale cloud shadow estimates axe combined with local scale radiometric 
measurements the result may be biased, since locally the cloud coverage and shadows can 
differ substantially. Clouds can block the direct sunlight, but the sky may be clear directly 
above, and vice versa. 

The cloud mapping capabilities of the VIL enable the determination of cloud shadows 
on the ground. Cloud shadows are restored by geometrically projecting the cloud returns 
onto the ground using an actual sun angle. The cloud shadow estimates can be calculated 
for areas as small as ~ 2000 m 2 . Local cloud shadow and cloud coverage estimates of the 
VIL can be used to determine the cloud effects on local radiometer measurements. Then, 
the VIL cloud coverage estimates can be used to estimate the area-averaged radiation. 
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3. Results and Discussion 


Figure 14 summarizes the CBL mean depths, cloud base altitudes, cloud top altitudes, 
and cloud coverages derived from the VIL boundary layer volume scans recorded during 
the 1989 FIFE. The results are filtered with 15-minute long filters: the CBL mean depths 
with a median filter, the cloud bases with a minimum filter, the cloud tops with a maximum 
filter, and the cloud coverages with an average filter. 

When clouds are connected to the CBL, the CBL mean depth estimate should be 
interpreted as the altitude of maximum cloud echo variation, as for example, on July 30 
from 9:00 to 10:00 CDT, when the CBL mean depth estimate varies widely between cloud 
base and top. In a case of multiple cloud layers, the cloud base altitude represents the base 
of the lowest cloud layer, while cloud top altitude represents the top of the highest cloud 
layer. That is the case on August 5 and August 11, when there is occasionally more than 
1 km difference between the cloud base and top altitudes. Clouds above ~4 km are not 
observed from the volume scans, since the maximum elevation angle of the boundary layer 
scans limits the altitude range. For the same reason, cloud coverage may be underestimated 
for cloud layers higher them 3.5 km. 

If cloud coverage is high, the cloud top altitudes may be too low due to missing lines 
of clear sight. On overcast days, such as August 1-3, the cloud top estimates are probably 
too low due to optically thick clouds blocking the backscatter profiles. For example, on 
August 2 from 10:00 to 12:00 CDT, when the cloud coverage is close to 100%, the cloud 
heights are limited by the backscatter signal penetration depth in clouds. When the cloud 
coverage is less than 70%, higher cloud heights were observed indicating that the method 
works successfully. 
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4. Conclusions 


Automatic methods for obtaining cloud coverages, cloud base altitude, cloud top alti- 
tudes, and convective boundary layer mean depths provide measurements which are aver- 
aged over 25-70 km 2 . All VIL data recorded during the 1989 FIFE program axe analyzed, 
providing information about boundary layer parameters in various atmospheric conditions. 
The accuracies of CBL mean depths and cloud base altitudes axe limited mainly by the 
vertical resolution of the VIL. 

The automatic CBL mean depth estimates correlate with the manually inspected 
estimates from RHI scans. Both the automatic and manual method for CBL mean depth 
estimates may have laxge uncertainties when the CBL is lower than 200 m, since the lidar 
frequently detects only the plume tops. The definition of the CBL mean depth becomes 
unclear in the presence of convective clouds. The automatic convective boundaxy layer 
mean depth estimate is valid if the coverage of convective clouds is less than ~10%. If 
more clouds are present, the maximum aerosol backscatter variance represents the altitude 
of maximum cloud echo variance, which usually is higher than the CBL mean depth. 

In this study, the maximum altitude of cloud detection is limited below ~ 3.5 km by 
the maximum elevation angle of the boundaxy layer volume scan. The root-mean-square 
difference between automatic and manual cloud base estimates is ~67 m. The cloud 
top estimates are usually too low in overcast conditions due to lidax signal extinction in 
optically thick clouds. When a 15-minute long maximum filter is applied to the cloud top 
estimates, there usually axe enough clear lines of sight through the cloud layer for reliable 
cloud top determination up to 70% cloud coverage. 
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Tables: 


Table 1 . VIL Specifications 

Transmitter 

Wavelength 

1064 nm (Nd:YAG) 

Average Power 

20 W 

Repetition Rate 

30 Hz 

Receiver 

Telescope Diameter 

50 cm 

Optical Bandwidth 

1 nm 

APD Quantum Efficiency 

~35 % 

Range Resolution 

15 m 

Maximum Angular Scan Rate 

20°s -1 

Data Rate 

~ 500 MB/h 

Data Acquisition 

Controlling Computer 

VAX 11/750 

Real-time Graphics 

Stellar GS1000 

Optical Disk 

C.Itoh OPCA-11 (2.6 GB disks) 

Data Processing 

IBM RS6000 Model 550 




Figure Captions: 


Fig. 1. The VIL scan pattern employed during the 1989 FIFE. A volume scan con- 
sists of an azimuthal sequence of elevation angle scans. The shaded areas illustrate 
aerosol structures. 

Fig. 2. An RHI scan, the variance of the filtered backscatter signal as a function 
of altitude (white curve), and the convective boundary layer mean depth estimate 
(black dashed line) on August 9, 1989, at 11:03 CDT. Aerosol structures are shown 
as brighter areas; cleaner air penetrating from free atmosphere is darker. The CBL 
top undulations are clearly visible. 

Fig. 3. Comparison between manually (solid circles) and automatically (open cir- 
cles) determined boundary layer mean depths on July 28, 1989. Open squares mark 
CBL mean depth estimates determined from radiosonde-based potential temper- 
ature profiles at 8:01 and 9:36 CDT. All RHI scans from 9:33 to 9:36 CDT were 
carefully reanalyzed to explain the difference between the VIL measurements and 
the radiosonde-based measurement. The solid circle at 9:33 CDT shows the average 
CBL depths obtained from individually inspected RHI scans; the triangle shows the 
average altitude of the deepest downward streams of clear air; the diamond shows 
the average altitude of highest plume tops. The vertical bars show minimas and 
maximas of the estimates. The plume top estimates indicate that the radiosonde 
profile at 9:36 CDT was measured through a high thermal plume top. 

Fig. 4. Potential temperature profiles measured with radiosondes on July 28, 1989, 
at 8:01 and 9:36 CDT. 



Fig. 5. Automatically vs. manually determined convective boundary layer mean 
depth with low cloud coverage from July 26 to August 11, 1989. Boxes a and c 
mark points where the CBL mean depth was hard to determine using the manual 
method due to low contrast between the boundary layer aerosols and free air; box 
b marks points which were manually misinterpreted. After reinterpretation, the 
manual estimates in these boxes are within 100 m of the automatic estimates. 

Fig. 6. An RHI scan image, a variance of the filtered backscatter signal as a function 
of altitude (white curve) and the CBL mean depth estimated from the variance 
profile (dashed black line) on August 2, 1989, at 12:00 CDT. Bright structures on 
the RHI scan represent cloud echoes. In this case, clouds make estimation of the 
CBL mean depth uncertain. 

Fig. 7. A VIL backscatter profile on July 26, 1989, at 10:48 CDT. Strong signal 
gradients characterize cloud echoes. The elevation angle of the profile is 15°. The 
irregularities below the first cloud return are due to boundary layer aerosols. The 
insert figure shows aerosol fluctuations expanded by a factor of ten. 

Fig. 8. An RHI scan, a histogram of the first cloud echoes (white line), and 
estimated cloud base (black dashed line) on August 2, 1989, at 12:00 CDT. Bright 
structures represent cloud echoes. Histogram peaks near the ground are caused by 
echoes from dust clouds and ground objects. 



Fig. 9. Cloud base altitude estimates using the manual method vs. automatic 
method from July 26 to August 11, 1989. Boxes a— d mark regions where the manual 
method missed the lower clouds due to sparse sampling. The highest sample in box 
c is due to an unusually high dust cloud, which caused the automatic method to 
fail. Box e marks a sample where aerosol swelling from a cloud was misinterpreted 
as a cloud in manual estimation. Box f marks a sample where the top of a thermal 
plume was misinterpreted as a cloud. 

Fig. 10. An RHI scan, a histogram of the last cloud returns (white line), and 
estimated cloud top (black dashed line) on August 2, 1989, at 12:00 CDT. 

Fig. 11. Cloud top altitude estimates using the manual method vs. automatic 
method from July 26 to August 11, 1989. Box a marks a point where a dust cloud 
at 500 m caused the automatic method to fail. Boxes b and e mark points where 
a human error caused spurious values in the manual estimates. Boxes c, d, and f 
mark points where a higher cloud layer was missed by the manual method, since 
the vertical range of the display was set too small. 

Fig. 12. Horizontal map of VIL cloud echoes on the ground on July 26, 1989, 
at 11:00 CDT. Grey and white areas indicate coverage by backscatter profiles and 
cloud echoes, respectively. 

Fig. 13. Comparison between a VIL cloud map (rectangle at upper-right corner) 
and an AVHRR visual band image (background) on July 28, 1989, at 15:00 CDT. 
The AVHRR image is enhanced to show clouds as white and the ground as black. 
The arc represents the VIL scanning area. 



Fig. 14. Summary of CBL mean depths (lines with open circles), cloud bases (lines 
with triangles up), cloud tops (lines with triangles down), and cloud coverages (solid 
line) for the 1989 FIFE. 
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Abstract 


This paper presents a study of the vertical profiles of horizontal mean winds calculated 
from data obtained with the University of Wisconsin Volume Imaging Lidar (VIL) during 
the 1989 FIFE program. Wind profiles are determined from cross correlation functions 
between subsequent Constant Altitude Plan Position Indicator planes calculated from vol- 
ume scans of the convective boundary layer. An objective technique to identify unreliable 
measurements is developed. An error analysis based on determining the errors due to 
random noise in data is performed. The VIL wind profiles are compared with radiosonde, 
aircraft-based, and ground-based measurements. Based on internal consistency tests, we 
estimate that in the convective boundary layer the hourly averaged VIL wind estimates 
which have passed the reliability analysis, have root-mean-square errors 0.03-0.1 ms -1 in 
speed and 0.4-2. 0° in direction. 
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1. Introduction 


Mesoscale wind profiles are needed in a wide variety of meteorological applications. 
Wind profiles are traditionally measured using methods which provide wind speed and di- 
rection at certain points. However, point measurements in the convective boundary layer 
do not generally provide representative area averaged wind estimates. The point mea- 
surements can be biased, for example, by large eddies, roll vortices, and topography. The 
averaging time cannot be extended to dampen temporal fluctuations without producing 
errors due to the changing value of the mean wind. Therefore, representative area average 
measurements are required. 

The three-dimensional aerosol detection capability of the University of Wisconsin 
Volume Imaging Lidax (VIL) [Eloranta and Forrest, 1992] along with correlation techniques 
provide vertical profiles of horizontal area average winds [Schols and Eloranta, 1992]. The 
volume scans of the VIL show the three-dimensional evolution of the aerosol distribution 
in a 15 km long, 40° wide, and 9-15° high sector consisting of ~5500 lidax backscatter 
profiles with a 15 meter range resolution. A three-minute scan time ensures that most of 
the aerosol structures appearing in a volume scan are coherent in the next scan. This paper 
studies the accuracy of wind profiles calculated from the data obtained with the VIL during 
the 1989 FIFE (First ISLSCP (International Satellite Land Surface Climatology Project) 
Field Experiment) program [Sellers et a 1., 1992] between July 26 and August 11. The 
original wind profiling algorithm introduced by Schols and Eloranta (1992) is refined. 

The wind profiling method is based on determining the movements of aerosol struc- 
tures from subsequent Constant Altitude Plan Position Indicator (CAPPI) scans of the 
VIL. The CAPPI scans represent two-dimensional horizontal maps of the convective bound- 
ary layer (CBL) aerosols at different altitudes. The aerosol movements axe determined by 
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calculating spatial two-dimensional cross correlation functions (CCF) between subsequent 
CAPPI scans at each altitude. The position of the CCF maximum peak indicates the 
mean wind speed and direction average over ~ 70 km 2 VIL scanning area. 

At the top of the convective boundary layer, where the backscatter signals have a good 
contrast between aerosol structures and clearer air, the wind estimation can be done di- 
rectly from the two-dimensional cross correlation function between two subsequent CAPPI 
scans. Above the top of the CBL and neax the surface, however, the situation is more 
complicated. Above the top, lidar signal extinction on clouds causes correlations between 
shadowed portions of the signal. Near the ground, aerosol plumes anchored to particular lo- 
cations cause stationary aerosol structures which lead to strong zero space lag correlations. 
Intense structures appearing only on one CAPPI may cause strong random correlations. 
These conditions prompt a need for preprocessing the backscatter signal and CAPPI scans 
before calculating the cross correlation functions. 

To eliminate the effects of variable atmospheric attenuation, scan angle dependent 
background signals, and shot to shot variations in energy normalization, the lidar backscat- 
ter profiles are median high pass filtered. Signals above optically thick clouds are not used 
to prevent correlations between cloud shadows. Since wind moves aerosol structures during 
scanning, the position of each lidar backscatter profile is shifted upwind to provide aerosol 
information at the time of the first shot in the volume scan. Previous hourly averaged 
results are used as a priori information and the wind profiling procedure is iterated if the 
subsequent results differ substantially. 

A temporal median of five subsequent CAPPI scans is subtracted from each CAPPI 
scan to prevent interference from stationary structures. The CAPPI scans are then his- 
togram equalized [Hall, 1974] to prevent any single intense aerosol structure from domi- 
nating the cross correlation function. In the earlier study [Scho/s and Eioranta, 1992], the 
histogram equalization was performed with 64 intensity classes. Because too few classes 


3 


will blur aerosol structures and widen cross correlation functions, the number of classes is 
limited only by the digital levels of the data in this study. Finally, the average intensi- 
ties axe subtracted from the CAPPI scans to reduce zero space lag correlations. Figure 1 
represents an example of a filtered CAPPI scan, which is ready for correlation calculations. 

A two-dimensional cross correlation function between two subsequent CAPPI scans 
at the same altitude is calculated using the Fast Fourier Transformation with zero padded 
data to prevent end effects [Press et a/., 1988]. The cross correlation function is divided 
by the covariance of the corresponding CAPPI scans providing a correctly weighted CCF 
for the time-average calculations. The wind estimate is determined from the maximum 
correlation peak in three steps: 

1) locating the global maximum correlation peak 

2) determining the local maximum of the global maximum peak 

3) calculating the wind speed and direction 

The global maximum correlation peak was determined from a spatially low pass fil- 
tered correlation function in the earlier study [Schols and Eloranta, 1992]. The optimal 
length for this filter, however, depends on the noise level and aerosol structure formations, 
which both vary on different atmospheric conditions. To reduce time and space dependent 
parameters in the analysis, this study uses the following approach to locate the maximum 
correlation peak. First, the global maximum correlation amplitude is found. All regions 
with amplitude greater than 1/e of the global maximum correlation amplitude are located. 
The mass m of each region is calculated by summing all correlation amplitudes inside the 
region. The region with greatest mass m max represents the global maximum correlation 
peak. The masses of the regions are also used in the reliability analysis (see section 2.). 

The local maximum of the global maximum correlation peak position is determined 
by least-squares- fitting a two-dimensional polynomial function to a 5x5 domain around 
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the highest point of the global maximum correlation peak. The maximum point of the 
fitted function is used as the maximum peak position gaining sub-pixel resolution. Sub- 
pixel resolution is achievable since aerosol structure formations along with approximate 
Gaussian turbulence lead to smooth-peaked cross correlation functions. The fitted function 
is expressed as 


/(x,y) = a 0 + a^x + a 2 y + a 3 x 2 + a 4 y 2 + a 5 xy 


( 1 ) 


where x and y denote coordinates in correlation plane and a<’s are model parameters. The 
model parameters (ao, aj , . . . , as) are solved from a linear equation group 
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where ( x;,y* ) are the coordinates of the fit point i and cc/(x,y) is the cross correlation 
function amplitude at point (x,y). The 25x6 kernel matrix is marked as G, the model 
parameter vector as a, and the vector of cross correlation function amplitudes as d. The 
least-squares solution of the parameter vector a is [Menke, 1984] 
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where G T denotes the transpose of matrix G and [ ] -1 the inverse of a square matrix. 

After determining the model parameters, the position of the maximum point (x max ,y max ) 
of function /(x, y) is determined by solving equations 
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leading to 


df(x,y) 

dx 

9f(x, y) 

dy 


(4) 

(5) 


2 fl 4 ( Z ] — 0.5(12 

Zmaz = 9 4 

05 — 4a 3 a 4 
2 a 3 a 2 — 0.50-1 

Vmaz = 7 

ag — 4a 3 a 4 

Figure 2 represents an example of the fit to the maximum CCF peak. The top of 
the peak is wide and the points near the local maximum have almost the same intensity. 
The maximum point of the fitted function interpolates between pixels gaining sub-pixel 
resolution. The polynomial fit is less sensitive to noise than the bicubic spline fit used in 
the earlier studies [ScAols a ad Eloranta, 1992]; the least-squares fit slightly filters the peak, 
while the bicubic spline may have noise induced oscillations near the maximum peak. 

The wind speed and direction are determined from the position of the maximum peak 
of the fitted polynomial representing mean movements of the aerosol structures between 
the CAPPI scans. The wind speed u and direction 6 are given by 


( 6 ) 
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where r is the time separation between subsequent volume scans. 


The most effective way to calculate time- averaged wind estimates is to average cross 
correlation functions together, since random correlations vanish in the averaging process 
providing a less noisy correlation function to the wind calculation. Even a weak aerosol 
correlation peak dominates the cross correlation function after an appropriate averaging 
time. Figure 3 compares the results of two different averaging methods. The first method 
averages subsequent correlation functions leading to a consistent wind profile. The second 
method calculates wind estimates from the cross correlation functions and then averages 
the wind estimates. Some of the estimates of the second method have large fluctuations 
from adjacent points due to a few uncertain results with large errors which were averaged 
together with reliable results. All results presented in this study are generated by averaging 
correlation functions. 

The main purpose of this study is to evaluate the accuracy and reliability of the 
wind measurements calculated from the data obtained with the VIL during the 1989 FIFE 
program. A method is developed for estimating the root-mean-square (RMS) errors in 
determining the wind speed and direction. If turbulent mixing changes aerosol structures 
between CAPPI scans substantially or the signal is dominated by noise, the strength of the 
aerosol structure correlations may decrease to the level of random noise correlations. This 
causes uncertain wind estimates with large errors. This study incorporates an objective 
method to determine the probability that the chosen correlation peak is caused by aerosol 
movements between two subsequent CAPPI scans instead of random signal fluctuations. 
This analysis is employed in the accuracy calculations to reject results which are likely to 
contain large errors. Reliable wind results are then compared with measurements made 
with an aircraft based instrument, radiosondes, and automatic weather stations. Even if 
these traditional measurements had no noise, they may not represent well area average 
winds due to bias by large scale eddies and roll vortices. Therefore, the method for esti- 
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mating the RMS errors of the wind results is tested comparing the error estimates with the 
internal consistency of the VIL results. The 1989 FIFE program produced a large number 
of VIL wind profiles which are entered into the FIFE database [Strebel et a/., 1990]. Only 
examples of them are presented here. 

2. Reliability Estimation 


The major problem in the wind profiling method is the possibility that the selected 
maximum correlation peak results from correlations of random noise instead of aerosol 
structures. Normally, a cross correlation function between two CAPPI scans is dominated 
by a single peak. Sometimes, however, random noise fluctuations along with weak aerosol 
structure coherence between subsequent CAPPI scans lead to random correlation peaks, 
which are stronger than the aerosol correlation peak. This prompts a need to identify the 
spurious results, which are likely to be calculated from random correlation peaks. 

In this study, an objective method is developed to distinguish between reliable and un- 
reliable wind estimates by indicating the probability that the chosen maximum correlation 
peak represents the mean aerosol movements. Statistical studies indicated the intensity 
of wind estimate fluctuations from the hourly mean values correlated with the CCF max- 
imum peak strength compared to all substantial peaks in the cross correlation function. 
Probability P mai that the global maximum correlation peak represents the mean wind is 
given by 


Pmaz — 


m r 


LL m. 


( 11 ) 


where m max is the mass of the global maximum correlation peak and m, is the mass 
of the correlation peak i. The masses are calculated by summing the intensities around 
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the peak with amplitude higher than 1/e of the global maximum correlation amplitude. 
Pmax indicates the probability that the strongest correlation function peak is caused by 
the mean wind. If there is only one correlation peak, then the P max is 1. If there are 
numerous similar correlation peaks, P max converges to zero. 


Figure 4 shows scatter plots for P max vs. wind speed and direction for quarter-hourly 
averaged results on July 27, 1989. Most of the results have P max = 1 and are within the 
natural wind deviation of that day; the results with lower P max are more widely spread. 
In general, P max > 0.5 provide mostly good results, while P max < 0.5 are poor. 


The fraction of reliable VIL wind estimates as a function of averaging time interval and 
Pmax is shown in figure 5. Averaging the cross correlation functions effectively increases 
the fraction of reliable results. Even a fifteen-minute averaging more than doubles the 
ratio of good measurements to all measurements. The recovery ratio of hourly averaged 
winds in the boundary layer is about 76%. Most of the measurements with P maz between 
0.9 and 1.0 are still consistent. 
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3. Error Estimation 


The errors of the wind speed and direction are estimated by calculating the change in 
the maximum correlation peak position (x mai ,j/ max ) when the CCF has some additional 
noise. In this section, we determine the error of the fit parameters and use it to estimate 
error in position of the maximum point of the fitted function. The expected root-mean- 
square error in wind speed and direction are derived from the error in the position of the 
maximum correlation peak. 

Assuming the noise in the CCF plane is uncorrelated from point to point and the 
variance of the noise is a 2 , the covariance matrix of the estimated model parameters a eat 
solved from equation (3) is given by 


cov a„t] = c 2 n 
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Since the histogram equalization is a highly nonlinear operation, we have no informa- 
tion about the spatial distribution of noise. We have assumed that the noise is uncorrelated 
with a constant variance to make analytical formulation of the noise effects possible; how- 
ever, we cannot determine the accuracy of this approximation. Although, after histogram 
equalization, the uniform distribution of the signal and noise superposition suggests that 
the noise in the CCF is probably also uniformly distributed. 


The noise variance a 2 of the CCF is estimated from a noise plane. The noise plane 
is made from the CCF by reducing the aerosol correlation portion of the CCF using the 
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following procedure. First, an average of three vertically adjacent CCF planes axe sub- 
tracted from the CCF. Then, the global maximum correlation peak is masked by rejecting 
all pixels inside the maximum correlation peaks 1/e-region. Finally, noise variance a n 
calculated from the masked CCF. 

The mean-squared error in determining the maximum peak position is given by 


and 
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where a 2 ’s are error variances of model parameters a,- estimated from the model parameter 
covariance matrix [cov a e3 t] as 
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Finally, the root- mean- square errors a u and <jq for wind speed and direction estimates, 
respectively, are determined from 
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The direction, error is expressed in radians. These equations represent root-mean-square 
errors of the mean wind estimate due to position error of the maximum correlation peak 
from a noisy CCF. The error estimates represent the lower bounds of the error, since 
the previous error analysis does not take into account possible errors due to systematic 
deformation of aerosol structures between scans. Furthermore, the error variance may 
be poorly estimated due to nonlinear effects in the histogram equalization of the CAPPI 
scans. 


4. Results and Discussion 


4.1. Reliability analysis 


This section summarizes the results of the reliability analysis performed on all VIL 
wind profiles from the 1989 FIFE program. Figure 6 shows probability P ma x of quarter- 
hourly averaged wind measurements as a function of time and altitude on July 26, 1989. 
The contours show the wind speed increasing from 5 ms -1 to 8.5 ms -1 over the period. 
Even using 0.2 ms -1 contour spacing, the wind profiles are smooth and small details are 
visible. The measurements have good internal consistency, except for some small distur- 
bances with a smaller P max . The results are most reliable in the convective boundary layer. 
Above the convective boundary layer there are sometimes not enough aerosol structures 
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for correlation analysis, leading to less reliability. Turbulent mixing in the late afternoon 
sometimes produces uniform aerosol distribution that leads to less reliable wind estimates, 
which also are present in the Figure 6. In the convective boundary layer, randomly located 
decreased reliabilities are generally due to the breaks in the data record. This happened 
on July 26 several times after 14:00 CDT. These breaks are generally due to system ad- 
justments or safety-shutoff of the laser output beam when an aircraft was flying in the field 
of view. 

4.2. Comparisons with flight-leg averaged and point measurements 


To verify consistency between the VIL wind profiles and traditional measurements, 
the VIL results are compared with: 

1) data from eight National Center for Atmospheric Studies Portable Automatic 
Mesonet stations (PAMS) [Sellers et al., 1992] that provide half-hour averaged 
ground-based measurements 

2) radiosonde profiles measured by a team from Cornell University directed by 
Prof W. Brutsaert [Sellers et al., 1992] 

3) aircraft-based Canadian National Aeronautical Establishments Twin Otter 
[ MacPherson , 1988] flight-leg averaged measurements 

All comparison data is retrieved from the FIFE database [Strebel et al., 1990]. Figure 
7 provides an overview of the site locations and flight patterns on July 27, 1989. All 
measurements are done within 20 km distance from each other. All weather stations axe 
closer than 10 km from the VIL site, and an average 20 m lower than the VIL site. 

Figure 8 compares hourly averaged VIL wind profiles, radiosonde profiles, aircraft- 
based measurements, and hourly averaged surface measurements on July 27, 1989, from 
11:00 to 12:00 CDT. The cloud base altitude is ~ 900 m; the cloud top altitude is ~ 1500 m; 
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cloud coverage is ~ 35%. The VIL wind profiles vary smoothly over the mixed layer and 
the cloud layer. The P ma x indicates less reliability above the cloud top where the aerosol 
distribution is weaker, although the method still gives some consistent results. As expected, 
wind speeds decrease slightly towards the ground showing the effects of friction. The wind 
shows little directional shear in the convective boundary layer. 

The surface measurements show averages of eight weather stations over one hour pe- 
riod. Surface friction causes the averaged ground winds to be slightly slower and backed 
compared to the low level VIL profiles. For safety reasons, the VIL did not scan at the 
ground level causing a vertical gap between the VIL and PAMS measurements. The wind 
measurements vary ~ 1 ms -1 and 20° between stations due to topographic differences. Ev- 
idently, averaging several ground measurements together still does not achieve the stability 
of area average wind measurements. 

The discrepancy between radiosonde profiles and area-averaged VIL profiles is ex- 
pected, since the radiosonde is sensitive to local wind variations due to turbulence. The 
radiosonde profiles suggest fairly strong directional wind shear at the top of the mixed 
layer and near to ground, while the VIL and aircraft measurements do not. Small av- 
eraging time interval of the radiosonde makes it sensitive to local wind variations. The 
area-averaged VIL wind profiles and flight-leg averaged aircraft wind measurements axe 
more representative in a mesoscale sense than the radiosonde wind profiles. 

The aircraft flight-path averaged wind measurements agree with the VIL profiles. 
The largest deviations between aircraft-based measurements and VIL profiles occur at the 
mixed layer top and near the ground. This is probably due to horizontal roll vortices, which 
axe detected as trends in aircraft wind calculations. These eddies produce mostly vertical 
air motions in the middle of the mixed layer, where the VIL and aircraft measurements 
coincide better. 


14 


Figure 9 shows the correlation between VIL results with P maz = 1 and simultane- 
ous aircraft-based measurements. A least-squares-fitted line through the data points has a 
slope of 0.98±0.2 and an offset of 0.17±0.12 ms -1 for wind speed, and a slope of 1.01±0.01 
and an offset of 1.7°±1.9° for direction. The root-mean-square deviations between results 
axe ~ 1.0 ms -1 in speed and ~ 12° in direction. The aircraft flights on July 28, August 
2, and August 3, were directed through the middle of the VIL scan pattern. Thus, the 
sampling of the aircraft and lidar were more closely matched. These results show smaller 
deviations from the mean value in the VIL observations than in adjacent aircraft measure- 
ments. This suggests that the VIL results are statistically more stable than the aircraft 
measurements, as a result of area averaging. 

4.3. Accuracy analysis 


Natural wind deviations in point and line-averaged measurements dominates the dif- 
ferences between these measurements and the VIL observations. Therefore, the point and 
line-averaged measurements are not used as the reference for accuracy analysis of the wind 
profiling method. Calculated error estimates of the VIL wind profiles are compared to the 
local consistency of the measurements. The local consistency is calculated for each wind 
measurement by subtracting an average of its vertically adjacent wind measurements. In 
linear wind shear conditions, the local consistency represents measurement inaccuracy. 
Data from July 26, 1989, is selected (see Figure 10) for accuracy calculations, since the 
wind shear is more linear than on any other day in the 1989 FIFE program. 

Figure 11 compares the local root-mean-square deviations of the results from the aver- 
ages of vertically adjacent results and calculated RMS error estimates for hourly averaged 
wind profiles on July 26, 1989. The calculated error estimates represent RMS errors in 
wind determination due to noise in the CCF, which is dominating error source after we 
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have rejected the wind estimates with P mai less than 1. Therefore, they axe generally 
smaller than the local deviations, which are due to both the natural wind deviation and 
method inaccuracy. In the convective boundary layer (below 1400 m) the RMS errors 
in hourly averaged wind estimates are less than 0.14 ms -1 in speed and 2.1° in direction. 
Part of the local variation is due to real wind fluctuations. Thus, it provides a conservative 
estimate of the wind fluctuations. 


The calculated profiles show minimum RMS errors at 800 m: the contrast between 
scattering from aerosol structures and mixed air is very good just below the CBL top, 
since the clear air penetrating down to the CBL is not yet mixed with the boundary 
layer aerosols. Closer to the ground, more intense mixing smears the aerosol structures 
decreasing the signal to noise ratio of the CCF. Above the convective boundary layer, 
clouds block the signal decreasing the signal to noise ratio. 


Figure 12 shows RMS-error estimates of wind measurements as functions of the av- 
eraging time and the normalized altitude from July 26 to July 31, 1989. Error estimates 
are linearly interpolated into altitudes normalized by the convective boundary layer mean 
depths and averaged over time. In the CBL, the error estimates are almost constant with 
altitude, but above the CBL, they grow almost exponentially. In the convective boundary 
layer, the RMS errors are in average less than 0.2 ms -1 in speed and 3.0° in direction. In 
the middle of the convective boundary layer, where aerosol structures have good contrast 
against the background, the RMS-errors are in order of 0.03-0.05 ms -1 in speed and 0.3- 
0.8° in direction. Above the CBL, the RMS-errors vary between 0. 1-3.0 ms -1 in speed and 
1.0-30° in direction, since the correlations get poorer due to aerosol structures diminishing 
with altitude and clouds blocking the signal. The RMS-errors decrease as the averaging 
time interval increases, since the averaging of the CCFs reduces random correlations. 
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5. Conclusions 


This study provides enhancements to the original wind profiling algorithm [Schols and 
Eloranta, 1992]. A method is developed to identify correlation peaks, which are likely to be 
caused by random correlations. The analysis performed for the VIL wind profiles from the 
1989 FIFE data indicates that by the most restrictive reliability definition (i.e. P m ai = l)» 
~76% of hourly averaged wind estimates in the convective boundary layer are reliable. 

The reliable wind profiles correlate well with traditional wind measurements made 
with an aircraft-based instrument, radiosondes, and weather stations. The comparisons 
do not indicate any systematic errors in the VIL wind profiles. The VIL wind profiles 
are more stable than the traditional measurements, since both time and area averaging is 
performed, while traditional measurements are either time-averaged or line-averaged. 

The deviation between VIL wind profiles and traditional measurements is dominated 
by the natural wind fluctuations in the traditional measurements. Therefore, the errors in 
the VIL profiles are estimated from internal consistency. This is performed by comparing 
the wind estimates with the averages of the vertically adjacent estimates. Data from linear 
wind shear conditions are used to minimize deviations due to natural non-lineax wind shear. 
The local root- mean- square errors of the hourly averaged wind profiles with P ma x = 1 are 
0.03-0.1 ms -1 in speed and 0.4-2° in direction. Since the local deviations are partly due 
to real wind fluctuations, they provide conservative error estimates. 

The error in the interpolating fit due to noise in the CCF provides the RMS-errors of 
the wind estimates. These error estimates are slightly under estimated for two main rea- 
sons: 1) possible systematic aerosol structure deformation between the scans may shear the 
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correlation peak, and 2) the noise variance of the correlation function is hard to estimate, 
since the histogram equalization of the CAPPI scans changes non-linearly the statistical 
behavior of the noise. However, the error estimates provide lower bounds of the RMS-error. 
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Figure Captions: 


Fig. 1. A filtered, wind corrected CAPPI scan on August 11, 1989, at 13:02 CDT 
at 1100 m altitude. The solid line represents the first shot angle of the scan; the 
dashed line represents the last shot angle of the scan. The arrow indicates the wind 
direction; its length represents the amount of correction for air parcel movements 
between the first and last shot. The bright areas represent increased backscatter 
from aerosols. The correction for air parcel movements during scanning shears the 
scanning area. 

Fig. 2. The cross correlation function maximum peak (wire frame) on July 
28, 1989, from 12:00 to 12:03 CDT at 500 meters altitude, fit points (diamonds), 
and the peak determined from the fit (rectangle). The horizontal axes represent the 
distance from the maximum point in pixels (one pixel equals 47 m). The correlation 
coefficient values are plotted on the vertical axis. 

Fig. 3. Wind profiles after different averaging methods on August 9, 1989, from 
10:30 to 10:45 CDT. The solid line of solid circles is a wind profile calculated from 
the average of five subsequent cross correlation functions. The dashed line with open 
circles is a wind profile calculated by averaging five individual wind profiles, which 
are estimated from subsequent cross correlation functions. The vertical dashed line 
at 500 m marks the convective boundary layer mean depth. 


Fig. 4. P max vs. speed (left) and direction (right) for quarter-hourly averaged 
wind estimates on July 27, 1989, from 8:00 to 15:00 CDT. All measurements with 
Pmax = 1 are within natural deviation of the wind; the most widely fluctuating 
wind results have P max < 0.5. 

Fig. 5. Fraction of reliable results as a function of averaging time interval and P max 
from July 26 to August 11, 1989. Open circles show fraction of results with P m ax = 1 
and altitudes ranging from 50 m to 2000 m; black rectangles show measurements in 
the convective boundary layer with P max = 1; open diamonds show measurements 
in the convective boundary layer with P max > 0.9. Averaging of CCF increases the 
fraction of reliable results by reducing random correlations. 

Fig. 6. Quarter-hourly VIL wind speeds (contours), mean boundary layer (solid 
line), and cloud base (dash-dotted line) plotted as functions of time and altitude on 
July 26, 1989. The shaded areas indicate P ma x is less than 1. In general, the results 
are reliable from the ground to 500 m above the cloud base. There are gaps in the 
data record at ~ 14:00, 16:00, and 17:00 CDT causing P max to decrease. Mixing 
of the boundary layer decreases the aerosol structure contrast in the volume scans 
causing occasional less reliable estimates in the middle of the mixed layer at 13:00, 
14:40, and 15:30 CDT. 



Fig. 7. Site locations and aircraft flight paths on July 27, 1989, from 11:00 to 
12:00 CDT. The 40° sector outlines the VIL scanning area; the intersection of 
the rectangle and the sector (bold polygon) represents the area where the VIL wind 
profiles are calculated. Lines show aircraft flight paths; PAMS locations axe marked 
with +; the radiosonde launch site is marked with *. The air parcels in the VIL 
scanning volume are measured by the aircraft about half hour later. The radiosonde 
launch site is about six kilometers north of the VIL site. The air parcels measured 
by the VIL missed the radiosonde site by about three kilometers. 

Fig. 8. Hourly wind profiles of the VIL compared with radiosonde profiles, Twin 
Otter aircraft-based measurements, and PAMS ground measurements on July 27, 
1989, from 11:00 to 12:00 CDT. The VIL results are marked with solid black sym- 
bols; error bars represent estimated RMS errors. Square symbols represent the 
radiosonde profiles which are vertically averaged to match the vertical spacing of 
the VIL results. The aircraft measurements are separated in East- West (open di- 
amonds) and South-North (solid diamonds) directed flights and connected in tem- 
poral order starting from the lowest measurement. Hourly averaged measurements 
from eight PAMS are averaged together (solid triangle); the error bars indicate the 
smallest and greatest individual PAMS measurements. 

Fig. 9. Comparison between VIL and aircraft-based wind speeds (left) and 
direction (right) measurements from the 1989 FIFE, from July 26 to August 11. 



Fig. 10. Hourly averaged wind profiles on July 26, 1989, between 11:00 and 17:00 
CDT. Only results with P max = 1 are plotted. The mixed layer mean depth rises 
rapidly to ~ 1000 m before 11:00 CDT when the clouds start forming. The base 
of convective clouds varies between 1000 m and 1800 m during the measurement 
session. The wind speed shear is linear except in the first 200-300 meters and at 
the top of the CBL. Error bars show estimated root-mean-square errors. 

Fig. 11. Local deviations of hourly averaged wind speed and direction (open 
circles) and calculated RMS error estimates (black circles) on July 26, 1989, between 
11:00 and 16:00 CDT. The error bars indicate the minimum and maximum error of 
hourly averaged wind estimates during 5-hour observation session. The calculated 
error estimates are probably too small, since only error in the interpolating fit is 
taken into account. On the other hand, the local deviations provide conservative 
error estimates, since they include also some real wind fluctuations. 

Fig. 12. Root-mean-square error estimate profiles of reliable results as functions of 
averaging time interval and normalized altitude from July 26 to July 31, 1989. The 
altitude is normalized with the convective boundary layer mean depth. Note the 
logarithmic horizontal axes. The profiles show RMS errors of wind measurements 
with P max = 1 averaged over a 5-hour period. Errors of 3-minute (open circles) and 
hourly average (black circles) profiles are shown. The RMS errors decrease as the 
averaging time interval increases. The minimum errors are obtained in the middle 
of the convective boundary layer. Below it, turbulent mixing smears the aerosol 
structures decreasing the signal to noise ratio of the CCF. Above it, the clouds 
block part of the data also decreasing the signal to noise ratio. 


-5.5 km 0 km 5.5 km 



Fig. 1 



0.18 

0.16 

0.14 

-4 



4 


Fig. 2. 


Altitude (km) 



Wind Speed (m/s) 


Fig. 3. 










Altitude (km) 



Time (CDT) 


Fig. 6. 



Altitu 



Speed (m/s) 



I . W IU 


7T 



0.5 

0.0 


Fig. 8. 








Altitude (km) Altitude (km) 




Altitude (km) 



Speed Deviation (m/s) Direction Deviation (deg.) 


Fig. 11. 


Altitude (km) 




Normalized Altitude (z/Zi) 



Fig. 12. 





Calculation of Divergence and Vertical Motion 
from Volume-Imaging Lidar Data 


P. W. Young 
Department of Physics 
University of Wisconsin - Platteville 
Platteville, WI 53818 


and 


E. W. Eloranta 

University of Wisconsin - Madison 
Deparment of Atmospheric and Oceanic Sciences 
1225 West Dayton Street 
Madison, WI 53706 



2 


ABSTRACT 

Boundary layer divergences are calculated from a sequence of images of the 
inhomogeneities in aerosol backscattering observed with Volume-Imaging Lidar. Two- 
dimensional spatial cross correlations are calculated between successive horizontal plane 
maps, 10 km x 5 km, which are interpolated from the lidar returns. The first map represents 
the actual aerosol pattern; the second map is generated with an assumed divergence to remove 
the effects of that divergence from the image. The divergence is determined from an 
optimization of the correlation versus this assumed divergence. The results from this process 
are consistent with predictions for divergence over the spatial scale measured, and the 
apparent resolution of the results agrees with the image resolution. 
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1. INTRODUCTION 

Horizontal divergence and vertical velocity are important parameters in the modeling of 
atmospheric processes at all scales. However, these quantities are very difficult to measure. 
On the synoptic scale vertical motion is important to weather prediction, climate modeling, 
and the study of general circulation. At present the most reliable measures of large-scale 
vertical motion come from analyzed wind fields produced by the National Meteorological 
Center and the European Center for Medium-Range Weather Forecasts [Lambert, 1989]. 
Radar wind profilers offer the potential for direct measurements of these large-scale motions, 
but their accuracy has not been fully established [Gage, 1990]. 

Mesoscale divergences and vertical motions are important in the study of convective 
storm systems and other boundary layer convergence processes. Single- : and multi- Doppler 
radar systems have been used to measure average divergence over scales of 15 to 80 km (700 
to 20,000 km 2 ) [Browning and Wexler, 1968; Rabin andZmic, 1980]. The divergence and 
other turbulence terms were determined from harmonic analysis of the azimuthal Doppler 
velocity field obtained with the scanning Doppler radar. The maximum fair-weather 
divergences measured over the above scales were of the order 10' 5 s' 1 . Vertical velocities 
were calculated from the divergences. 

Less effort has been put into measuring smaller scale divergences and vertical velocities, 
although they influence the larger scale motions and are very important in boundary layer 
processes. Also, these smaller scales encompass subjects such as flow over obstacles and 
pollution studies. Doppler lidar data have been used to determine divergence and other 
turbulence parameters over scales between 150 m and 12 km [ Eberhard, et al, 1989; 
Gultepe, et al, 1990; Gal-Chen, et al, 1992]. The divergences were calculated using the 
VAD technique developed for Doppler radar. The turbulent fluxes were computed from the 
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variance of the radial velocities. Small-scale vertical motions in the boundary layer have also 
been measured directly with radar wind profilers [Ecklund, et al 1988] and tower-mounted 
anemometers [Skibin, et al, 1985]. The measurements were averaged over short time 
intervals to provide the small-scale motions. Divergences at ground level have been 
measured over scales from 300 m to 6 km (0.04 to 15 km 2 ) using an optical propagation 
technique [Kjelaas and Ochs, 1974; Fritz and Wang, 1979]. The divergence was calculated 
from transverse velocities measured from the optical propagation along three legs of a 
triangle. The maximum values of the divergence measured over these areas ranged from 10' 2 
s' 1 for the smaller area down to 10' 3 s' 1 for the larger area. Figure 1 summarizes these clear- 
weather divergence measurements over the micro- and meso-scales. 

This paper describes a method for determining the spatially averaged horizontal 
divergence in the boundary layer from Volume-Imaging Lidar (VIL) data. The VIL produces 
a time sequence of images of the aerosol structure in 30° to 40° sectors of the boundary layer 
out to 15 km (60 km 2 ). Sequential images are correlated to determine the vertical profile of 
the horizontal divergence over the area of the images. From Figure 1 the maximum fair- 
weather divergences over the VIL area are expected to be of the order 10" 4 s' 1 . 

2. METHOD OF DATA ANALYSIS 

The time evolution of the atmospheric aerosol distribution in a particular volume of the 
boundary layer is due to the wind and turbulence, unless there are sources or sinks inside the 
area. The mean wind translates the distribution. Spatial variations in the wind, including the 
divergence, alter the pattern of the inhomogeneities in the distribution. VIL produces a time 

sequence of maps of this distribution from the aerosol backscattering of the laser beam. For 
wind and divergence calculations the data are processed to produce a vertical stack of 
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horizontal images which represent the aerosol distribution at constant altitude. The evenly 
spaced grid points in the images are produced by linear interpolation of the lidar returns. 
Each data point in the rectangular grid is the weighted average of all the lidar returns located 
in the eight cells surrounding the point. The weighting factor is inversely proportional to the 
distance from the grid point to the coordinates of the lidar return. The interpolation process 
also includes a first order correction for the image distortion caused by the mean wind over 
the scan time [Schols and Eloranta, 1992]. 

Calculations of vertical profiles of the mean horizontal wind from cross correlations of 
VIL images have been performed by Schols and Eloranta [1992] and Piironen and Eloranta 
[this issue]. The two-dimensional spatial cross correlation function (CCF) is calculated 
between two successive images in the same horizontal plane. The coordinates of the 
maximum of the CCF represent the mean displacement of the inhomogeneities in the 
distribution between the two images. This displacement divided by the time interval between 
the two images gives the mean horizontal velocity over the area. Time averaging of the mean 
velocity is accomplished by averaging the CCF for several time intervals before computing 
the velocity. These calculations are performed on image planes at different altitudes to 
construct a vertical profile of the horizontal velocity. 

This method for calculating the mean wind over an area averages out the effects of the 
spatial variations in the wind. However, the horizontal divergence can be determined using 
the CCF method by incorporating the effects of the divergence on the two-dimensional 
distribution before calculating the CCF. The divergence in the horizontal wind field 
stretches (or compresses) the inhomogeneities in the aerosol distribution over time. The 
aerosol in a rectangle of area A = L x x L y at time tj is translated by the mean wind and 
spread by the divergence into the following area at time 12 = 1!+ At: [Hess, 1959] 



6 


A' 





(1) 


where 


x,y = position coordinates in area, with x along the sector axis and 
y transverse to the axis; 

u,v = velocity components along x and y, respectively; 

L x ,L y = image dimensions at time t t ; 

At = time between images. 

Therefore, in the VEL image at time t 2 the aerosol patterns are not only shifted by the mean 
wind, they are also spread out (or compressed) horizontally by the divergence. The effect 
of the divergence can be removed from the two-dimensional VIL image at time t 2 by 
recalculating an image of size A from the VIL return data from area A’ centered on the mean 
wind displacement. The interpolation process is repeated with a grid size based on A’ while 
the grid size in the final image is left consistent with A. This effectively compresses (or 
stretches) the inhomogeneities in the image back into their original shape. The maximum of 
the CCF between this "corrected" image at t 2 and the actual image at tj should be greater 
than the maximum of the corresponding CCF between the actual images at both times. Thus, 
the value of the horizontal divergence in the time interval between two VIL images can be 
determined by optimizing the maximum of the CCF between the two images versus du/dx and 
dv/d y. 

This process has been applied to lidar data taken with the University of Wisconsin VIL 
[Eloranta and Forrest , 1992]. The data were taken in a 40° sector out to a range of 15 km 
and up to an altitude of 2 km. For the divergence calculations the original image area was 
a 10 km x 5 km rectangle between the ranges of 4 km and 14 km. The images of this area 
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were processed into 200 x 100 data arrays. The vertical spacing between image planes was 
50 m. The time interval between successive images was 180 to 200 s. 

The maxima of the CCFs between successive images (corresponding to the mean wind 
displacement) ranged from about 0. 1 to 0.4 for the data analyzed. The image at time t 2 was 
taken from the same area as the image at tj so the CCF maximum was displaced from the 
origin by the mean wind. The value of the CCF maximum was determined from a two- 
dimensional, 2nd order polynomial fit to the CCF data about the peak. The correlation peak 
is not a 2nd order polynomial, but the polynomial is an acceptable approximation over a small 
square of grid points around the peak. When the actual image at time t 2 was replaced by an 
image processed from the area adjusted for the divergence according to equation (1), the CCF 
maximum improved as seen in Figure 2. The amount of improvement ranged from 0% 
(when the divergence was zero) to about 10%. 

Optimization of the CCF maximum was carried out by calculating the CCFs between the 
actual image at tj and images at t 2 adjusted for an array of assumed 3u/3x and 3v/3y in 
equation (1). The increments in 3u/3x and 3v/3y were chosen to match the image resolution, 
50 m. The images were altered symmetrically so the incremental changes in L x and L y were 
100 m, corresponding to increments of 5 x 10' 5 s' 1 and 10 x 10‘ 5 s' 1 for 3u/3x and dv/dy, 
respectively. Figure 3 shows a plot of these calculations. The optimum values for 3u/3x and 
3v/3y were determined by fitting a two-dimensional, 2nd-order polynomial in 3u/3x and 
3v/3y to the graph. The divergence was identified as the sum of these optimum values. 
Time averaging of the divergence was achieved by performing these optimization calculations 
on time-averaged CCFs (i.e. the average of several sequential CCFs). 

This CCF technique for calculating divergence is limited by the resolution in the two- 
dimensional images, 50 m x 50 m. Stretching the image by this amount produces 3u/3x and 
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dv/dy values of 2.5 x 10* 5 s' 1 and 5 x 10* 5 s' 1 , respectively. These values provide a rough 
estimate of the limiting resolution for the technique. The resulting resolution in the 
divergence is ~6 x 10' 5 s' 1 . This limit is considerably larger than the limitations on mean 
wind calculations from CCFs discussed in Piironen and Eloranta [this issue]. They report 
mean wind uncertainties of less than 0. 1 m/s using CCF calculations over approximately the 
same area. This uncertainty in the mean wind corresponds to a divergence of less than 
2 x 10' 5 s' 1 . 

The continuity equation was used to relate the gradient of the vertical motion to the 
horizontal divergence, 

(du + dv^j + dw _ 1 dp (2) 

dy) dz p dt 

where 

z = vertical coordinate; 
w = vertical component of velocity; 
p = density of fluid. 

For an incompressible fluid the gradient of the vertical wind is the negative of the horizontal 
divergence. Therefore, the average vertical motion over the area imaged by the VIL can be 
determined from the integral of the vertical profile of the divergences. 

3. RESULTS 

This section presents results of divergence and vertical motion calculations on VIL data 
taken as part of the First ISLSCP Field Experiment (FIFE) 1989 on the Konza Prairie 
research site near Manhattan, Kansas. The VIL scans were centered on azimuth 236°. Two 
data sets representing different wind conditions were chosen to test the divergence algorithm. 
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The first set was collected between 1145 and 1240 CDT on August 3, 1989. The sky was 
overcast with a stratocumulus clouddeck which was just beginning to break up during the data 
period. The lowest clouds extended down to 400 m at the beginning of period. By 1240 
CDT the lowest clouds had risen to 550 m. These low clouds interfered somewhat with the 
divergence analysis between 400 and 600 m. The wind was relatively constant beteen 10 and 
11 m/s at 220°. This wind was approaching the lidar 26° off the center angle of the scan. 
The wind profile midway through the time interval is shown in Figure 4. 

Because of the relatively high wind speed, the divergence calculations were performed 
on CCFs averaged over about 15 minutes (5 scans). Figure 5 shows the divergence profiles 
below the clouddeck at three times selected so the data averages do not overlap. Between 
about 150 m and 500 m the horizontal divergence appears relatively constant at an average 
of -6 x 10" 5 s' 1 with a standard deviation of 7 x 10' 5 s' 1 for the 24 data points. If the data 
at 150 m are omitted, the standard deviation is only 5 x 10' 5 s' 1 . Integrating the average 
divergence up to 500 m gives a 3 cm/s upward velocity at 500 m. The magnitude of the 
divergence is consistent with the expected divergence over 50 km 2 from Figure 1, and the 
standard deviation is approximately equal to the expected resolution, ~6 x 10' 5 s' 1 . 

Above 500 m and below 150 m the results fluctuate much more than the estimated 
resolution. The fluctuations above 500 m can be attributed to interference from the low- 
hanging clouds. The large fluctuations at low altitudes are consistent with the observed 
deterioration of low altitude velocity calculations from CCFs [Piironen and Eloranta, this 
issue]. 

The second data set was collected between 0817 and 1045 CDT on July 28, 1989. This 
day was mostly clear with some cirrus. A few scattered boundary layer clouds developed just 
after this period. The mixed layer increased from about 100 m at the start to about 750 m 
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at the end. The winds were weaker than on August 3, starting at about 3.5 m/s at 210° and 
increasing steadily to about 4.5 m/s at 195° by the end. The wind was approaching the lidar 
between 26° and 41° from the center angle of the scan. Figure 6 shows the wind profiles 
during the data interval. 

The divergence calculations for this day were performed on CCFs averaged over just 10 
minutes (3 scans). Figure 7 shows divergence profiles for three times over 30 minutes, the 
same time span observed on August 3. In contrast to Figure 5, the individual profiles are not 
constant below the mixed layer, and on this day the divergence changed significantly with 
time. At the lower altitudes the profiles shift from convergent to divergent within the half 
hour period. Over the entire 2 x h - -hour period, the divergence below the mixed layer varied 
from about -25 x 10' 5 to +25 x 10' 5 s' 1 . Figure 8 shows these variations at a single altitude. 
The divergence changed continuously for the full data period. 

Because the divergence was changing with both time and altitude, the standard deviation 
of all the divergence values would be greater than the resolution. Instead, the resolution is 
estimated from the standard deviation of the difference between divergences at successive, 
independent times which is 11 x 10' 5 s' 1 . These fluctuations in the difference correspond to 
an uncertainty of less than 8 x 10' 5 s' 1 on the divergence, comparable to the estimated 
resolution based on the image resolution. Furthermore, a temporal correlation calculation on 
the data at all elevations indicates a correlation time of 20 to 30 minutes for the divergence 
on this day. ^ 

However, it is also recognized that these observed correlations in the divergence might 
be produced by similarities in the aerosol maps which persist for 20 to 30 minutes. All the 
images involved in the divergence calculations represent the same data area, so the images 
are correlated. It is possible that this correlation in the images is producing the observed 
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correlation in the divergence. Under this assumption, the temporal variations in the 
divergence would merely represent the random variations due to the uncertainties in the 
measurements. The standard deviation of the data would be the resolution on the divergence. 
This value would represent the upper limit on the resolution for this technique. Figure 9 
represents this interpretation of the results for July 28. This figure contains plots of all the 
gradient ( du/dx and dv/dy) and divergence values measured below the mean mixed layer. 
The solid line on the divergence graph represents the divergence calculated from the CCF 
averaged over the entire 2 ! /2-hour time period. The data points in each graph all average out 
to about zero, and the scatter could be viewed as just the uncertainty in the measurement. 
The standard deviations for du/dx, dv/dy, and divergence are 7 x 10' 4 5 , 10 x 10' 5 , and 
11 x 10' 5 s* 1 , respectively. The scatter in dv/dy is greater than that in du/dx because the 
transverse dimension is smaller. The 11 x 10' 5 s* 1 standard deviation in the divergence 
represents the upper limit on the resolution for this technique. This value is about twice the 
resolution estimated from the image resolution, but it is still the same order of magnitude as 
the expected maximum divergence over the area. 


4. CONCLUSIONS 

Horizontal divergence values, spatially averaged over 50 km 2 , have been calculated from 
VIL data taken during FIFE. On an overcast day with medium winds the divergence below 
the clouds was relatively constant around -6 x 10' 5 s' 1 for over a half hour. On a clear day 
with light winds the divergence changed continuously between -25 x 10" 5 and +25 x 10' 5 s* 1 
with a correlation time around 20 minutes. These measured divergences are consistent with 
the expected scale variations presented in Figure 1 . The internal consistency of the results 
indicate that the uncertainty on the divergence is about 7 x 10‘ 5 s' 1 , a value consistent with 
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the spatial resolution of the images. This indicates that these results are as reliable as can be 
expected from the lidar data available. 

Even if this internal consistency of the results is rejected, the most conservative 
interpretation of the data yields an upper limit on the resolution for the technique of 1 1 x 10' 5 
s' 1 which is the standard deviation of all the results below the mean mixed layer. This upper 
limit is less than twice the resolution based on the image resolution, and it is still the same 
order of magnitude as the expected maximum divergence over a 50 km 2 area. With 
reasonable improvements in the data collection procedures, this resolution can be reduced 
enough to measure divergences on this scale. 

Further development of this technique for measuring divergence with VIL will require 
improvements in the VIL to improve the resolution. The present resolution is limited by the 
image size and resolution. Increasing the repetition rate of the VIL would improve both these 
limations. The rate could be increased from 30 Hz to 1 kHz without interference from wrap- 
around returns from distant clouds. About 300 Hz would allow for 360° scans with smaller 
azimuthal steps in the same time interval. With a 360° image out to 15 km the image grid 
could be increased to 20 x 20 km. The additional factor of three in the repetition rate would 
also allow for smaller azimuthal steps, thus decreasing the image resolution to about 20 m. 
These two factors combined would provide a divergence resolution of about 5 x 10" 6 s' 1 . The 
increased image area and resolution would also improve the image correlation, possibly 
improving the resolution even more. With this level of performance the VIL could be used 
to monitor even weak divergences over the 400 km 2 area, thus complementing the existing 
data on horizontal divergence. This area would also overlap with Doppler radar’s capabilities 
allowing for comparisons. 
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FIGURE CAPTIONS 


Fig. 1. Scale dependence of the fair-weather divergence measurements from the literature. 
The optical propagation data are plotted as the open triangle [Kjelaas and Ochs, 1974] and 
the solid triangles [Fritz and Wang, 1979]. Doppler radar measurements are plotted as solid 
squares [Rabin and Zmic, 1980] and open squares [Rabin and Zawadzki, 1984]. The solid 
circle represents lidar measurements by Boers, Eloranta, and Coulter [1984]. The open 
circles are the median values from an experiment involving several sensors reported by Stull 
[1988, after Vachalek, 1987]. The VIL measurements reported in this paper are spatially 
averaged over 50 km 2 . 

Fig. 2. A comparison of the CCF calculated between the actual images at both t t and t 2 (solid 
line) and the CCF calculated with a divergence-corrected image at time t 2 (dotted line). The 
largest cross correlation was found with du/dx = 10 x 10" 5 s' 1 and dv/dy = -10.0 x 10' 5 s' 1 . 
The graphs show the cross sections of the two-dimensional cross correlation function around 
the peak. The CCFs were averaged over 3 images (10 min.). 

Fig. 3. CCF maxima as a function of assumed divergence for the same lidar data represented 
in Figure 2. The divergence was determined by fitting a 2nd order polynomial to this graph: 
du/dx = 12.4 x 10' 5 s' 1 , dv/dy = -6.6 x 10" 5 s' 1 , divergence = 5.8 x 10' 5 s' 1 . 
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Fig. 4. Wind velocity profile at 1213 CDT on August 3, 1989. The velocity and direction 
were calculated from the shift in the location of the CCF peak. The CCFs were averaged 
over 5 images (17 min.). 

Fig. 5. Divergence profiles for 1157 CDT (circles), 1213 CDT (squares), and 1228 CDT 
(triangles) on August 3, 1989. The data from 150 m to 500 m appear constant at an average 
of -6 x 10' 5 s’ 1 with a standard deviation of 7 x 10* 5 s' 1 . 

Fig. 6. Wind velocity profiles at hour intervals on July 28, 1989. The velocity and direction 
were calculated from the shift in the location of the CCF peak. The CCFs were averaged 
over 3 images (10 min.). 

Fig. 7. Divergence profiles for 0922 CDT (circles), 0939 CDT (triangles), and 0954 CDT 
(squares) on July 28, 1989. The solid symbols with solid lines represent data below the mean 
mixed layer; the open symbols with dotted lines show additional data above the mixed layer. 

Fig. 8 Time evolution of the divergence at 200 m on July 28, 1989. Each point represents 
a 10 minute average; the points are spaced so the averages do not overlap. 
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Fig. 9. Scatter plots of all the gradient and divergence results below the mean mixed layer 
on July 28, 1989. The divergence data are the sum of the corresponding data for the 
gradients. The points were selected from profiles separated by 10 minutes so they represent 
independent CCFs. The solid line in the divergence profile shows the divergence values 
calculated from the CCFs averaged over the entire 2 x h hour period. 
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